# "Media's Influence on LGBTQ Support Across Africa" by Stephen Winkler

# Code to generate:
#   - A.5a-e (first differences of freedom house interaction by country)
#   - A.6a-e (first differences of KOF interaction by country)

rm(list=ls()) # clear environment
source("Rprofile.R") # setup Rprofile (load packages, etc.)

# Load data
myData <- readRDS("data/afrobarometer.rds")

# Transform variables:
myData$sexuality2 <- as.numeric(myData$sexuality2)
myData$tolerance_noLGBT2 <- as.numeric(myData$tolerance_noLGBT2)
myData$radio <- as.numeric(myData$radio)
myData$tv <- as.numeric(myData$tv)
myData$newspaper <- as.numeric(myData$newspaper)
myData$internet <- as.numeric(myData$internet)
myData$social_media <- as.numeric(myData$social_media)
myData$media_aggregate <- as.numeric(myData$media_aggregate)
myData$media_noradio <- as.numeric(myData$media_noradio)
myData$media_notv <- as.numeric(myData$media_notv)
myData$media_nopaper <- as.numeric(myData$media_nopaper)
myData$media_nointernet <- as.numeric(myData$media_nointernet)
myData$media_nosocialmedia <- as.numeric(myData$media_nosocialmedia)
myData$female <- as.numeric(myData$female)
myData$education <- as.numeric(myData$education)
myData$religiosity <- as.numeric(myData$religiosity)
myData$age <- as.character(myData$age) #age needs to go through character first
myData$age <- as.numeric(myData$age) #then to numeric to maintain real values
myData$income_water <- as.numeric(myData$income_water)
myData$urban <- as.numeric(myData$urban)

# Country Datasets
benin_data <- myData %>% dplyr::filter(country == "Benin")
botswana_data <- myData %>% dplyr::filter(country == "Botswana")
burkina_data <- myData %>% dplyr::filter(country == "Burkina Faso")
burundi_data <- myData %>% dplyr::filter(country == "Burundi")
cameroon_data <- myData %>% dplyr::filter(country == "Cameroon")
capeverde_data <- myData %>% dplyr::filter(country == "Cape Verde")
cdi_data <- myData %>% dplyr::filter(country == "Cote d'Ivoire")
gabon_data <- myData %>% dplyr::filter(country == "Gabon")
ghana_data <- myData %>% dplyr::filter(country == "Ghana")
guinea_data <- myData %>% dplyr::filter(country == "Guinea")
kenya_data <- myData %>% dplyr::filter(country == "Kenya")
lesotho_data <- myData %>% dplyr::filter(country == "Lesotho")
liberia_data <- myData %>% dplyr::filter(country == "Liberia")
madag_data <- myData %>% dplyr::filter(country == "Madagascar")
malawi_data <- myData %>% dplyr::filter(country == "Malawi")
mali_data <- myData %>% dplyr::filter(country == "Mali")
maurit_data <- myData %>% dplyr::filter(country == "Mauritius")
morocco_data <- myData %>% dplyr::filter(country == "Morocco")
mozamb_data <- myData %>% dplyr::filter(country == "Mozambique")
namibia_data <- myData %>% dplyr::filter(country == "Namibia")
niger_data <- myData %>% dplyr::filter(country == "Niger")
nigeria_data <- myData %>% dplyr::filter(country == "Nigeria")
stp_data <- myData %>% dplyr::filter(country == "São Tomé and Príncipe")
senegal_data <- myData %>% dplyr::filter(country == "Senegal")
sierra_data <- myData %>% dplyr::filter(country == "Sierra Leone")
southaf_data <- myData %>% dplyr::filter(country == "South Africa")
swaz_data <- myData %>% dplyr::filter(country == "Swaziland")
tanz_data <- myData %>% dplyr::filter(country == "Tanzania")
togo_data <- myData %>% dplyr::filter(country == "Togo")
tunisia_data <- myData %>% dplyr::filter(country == "Tunisia")
uganda_data <- myData %>% dplyr::filter(country == "Uganda")
zambia_data <- myData %>% dplyr::filter(country == "Zambia")
zimbab_data <- myData %>% dplyr::filter(country == "Zimbabwe")

sims <- 10000 # will use this to create FD

# Models

##################
# Media Aggregate
# VI by district
##################

# Don't do aggregate for now. The FD shift from 1 to 5 is not relevant for this media type. 

#################
# Radio
# VI by district
#################

radio.country <- #use this model to estimate result
  sexuality2 ~
  radio +
  media_noradio + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 | district)  #add VI for District

radio.country.nodist <- #use this model to create scenarios (dropping district)
  sexuality2 ~
  radio +
  media_noradio + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban

#Benin
mdata <- extractdata(radio.country, benin_data, na.rm = TRUE) #only extract data used in model
benin.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                                   nAGQ = 0 ) # optimizer for faster but less precise fit
benin.radio.pe <- fixef(benin.radio.result) # extract pe
benin.radio.vc <- vcov(benin.radio.result) # extract vcov
benin.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
benin.radio.scen <- cfName(benin.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
benin.radio.scen <- cfChange(benin.radio.scen, "radio", 
                          x = 5,
                          xpre=1,
                          scen=1)
benin.radio.simbetas <- mvrnorm(sims, benin.radio.pe, benin.radio.vc) #create simbetas
benin.radio.qoi <- logitsimfd(benin.radio.scen, benin.radio.simbetas, ci=0.95) #and get quantities of interest
benin.radio.qoi.df <- as.data.frame(benin.radio.qoi) %>% mutate(medium = 'radio', country = "benin")

#Botswana
mdata <- extractdata(radio.country.nodist, botswana_data, na.rm = TRUE) #use nodist model for botswana b/c there are no districts (i.e. no random effects)
botswana.radio.result <- glm(radio.country.nodist, family=binomial, data=mdata) #estimate the model
botswana.radio.pe <-botswana.radio.result$coefficients  # extract pe
botswana.radio.vc <- vcov(botswana.radio.result) # extract vcov
botswana.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
botswana.radio.scen <- cfName(botswana.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
botswana.radio.scen <- cfChange(botswana.radio.scen, "radio", 
                             x = 5,
                             xpre=1,
                             scen=1)
botswana.radio.simbetas <- mvrnorm(sims, botswana.radio.pe, botswana.radio.vc) #create simbetas
botswana.radio.qoi <- logitsimfd(botswana.radio.scen, botswana.radio.simbetas, ci=0.95) #and get quantities of interest
botswana.radio.qoi.df <- as.data.frame(botswana.radio.qoi) %>% mutate(medium = 'radio', country = "botswana")

#Burkina Faso
mdata <- extractdata(radio.country, burkina_data, na.rm = TRUE) #only extract data used in model
burkina.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                            nAGQ = 0 ) # optimizer for faster but less precise fit
burkina.radio.pe <- fixef(burkina.radio.result) # extract pe
burkina.radio.vc <- vcov(burkina.radio.result) # extract vcov
burkina.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burkina.radio.scen <- cfName(burkina.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
burkina.radio.scen <- cfChange(burkina.radio.scen, "radio", 
                             x = 5,
                             xpre=1,
                             scen=1)
burkina.radio.simbetas <- mvrnorm(sims, burkina.radio.pe, burkina.radio.vc) #create simbetas
burkina.radio.qoi <- logitsimfd(burkina.radio.scen, burkina.radio.simbetas, ci=0.95) #and get quantities of interest
burkina.radio.qoi.df <- as.data.frame(burkina.radio.qoi) %>% mutate(medium = 'radio', country = "burkina faso")

#burundi
mdata <- extractdata(radio.country, burundi_data, na.rm = TRUE) #only extract data used in model
burundi.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
burundi.radio.pe <- fixef(burundi.radio.result) # extract pe
burundi.radio.vc <- vcov(burundi.radio.result) # extract vcov
burundi.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burundi.radio.scen <- cfName(burundi.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
burundi.radio.scen <- cfChange(burundi.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
burundi.radio.simbetas <- mvrnorm(sims, burundi.radio.pe, burundi.radio.vc) #create simbetas
burundi.radio.qoi <- logitsimfd(burundi.radio.scen, burundi.radio.simbetas, ci=0.95) #and get quantities of interest
burundi.radio.qoi.df <- as.data.frame(burundi.radio.qoi) %>% mutate(medium = 'radio', country = "burundi")

#Cameroon
mdata <- extractdata(radio.country, cameroon_data, na.rm = TRUE) #only extract data used in model
cameroon.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
cameroon.radio.pe <- fixef(cameroon.radio.result) # extract pe
cameroon.radio.vc <- vcov(cameroon.radio.result) # extract vcov
cameroon.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cameroon.radio.scen <- cfName(cameroon.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
cameroon.radio.scen <- cfChange(cameroon.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
cameroon.radio.simbetas <- mvrnorm(sims, cameroon.radio.pe, cameroon.radio.vc) #create simbetas
cameroon.radio.qoi <- logitsimfd(cameroon.radio.scen, cameroon.radio.simbetas, ci=0.95) #and get quantities of interest
cameroon.radio.qoi.df <- as.data.frame(cameroon.radio.qoi) %>% mutate(medium = 'radio', country = "cameroon")

#Cape Verde
mdata <- extractdata(radio.country, capeverde_data, na.rm = TRUE) #only extract data used in model
capeverde.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
capeverde.radio.pe <- fixef(capeverde.radio.result) # extract pe
capeverde.radio.vc <- vcov(capeverde.radio.result) # extract vcov
capeverde.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
capeverde.radio.scen <- cfName(capeverde.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
capeverde.radio.scen <- cfChange(capeverde.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
capeverde.radio.simbetas <- mvrnorm(sims, capeverde.radio.pe, capeverde.radio.vc) #create simbetas
capeverde.radio.qoi <- logitsimfd(capeverde.radio.scen, capeverde.radio.simbetas, ci=0.95) #and get quantities of interest
capeverde.radio.qoi.df <- as.data.frame(capeverde.radio.qoi) %>% mutate(medium = 'radio', country = "cape verde")

#cote d'ivoire
mdata <- extractdata(radio.country, cdi_data, na.rm = TRUE) #only extract data used in model
cdi.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
cdi.radio.pe <- fixef(cdi.radio.result) # extract pe
cdi.radio.vc <- vcov(cdi.radio.result) # extract vcov
cdi.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cdi.radio.scen <- cfName(cdi.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
cdi.radio.scen <- cfChange(cdi.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
cdi.radio.simbetas <- mvrnorm(sims, cdi.radio.pe, cdi.radio.vc) #create simbetas
cdi.radio.qoi <- logitsimfd(cdi.radio.scen, cdi.radio.simbetas, ci=0.95) #and get quantities of interest
cdi.radio.qoi.df <- as.data.frame(cdi.radio.qoi) %>% mutate(medium = 'radio', country = "cote d'ivoire")

#gabon
mdata <- extractdata(radio.country, gabon_data, na.rm = TRUE) #only extract data used in model
gabon.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
gabon.radio.pe <- fixef(gabon.radio.result) # extract pe
gabon.radio.vc <- vcov(gabon.radio.result) # extract vcov
gabon.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
gabon.radio.scen <- cfName(gabon.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
gabon.radio.scen <- cfChange(gabon.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
gabon.radio.simbetas <- mvrnorm(sims, gabon.radio.pe, gabon.radio.vc) #create simbetas
gabon.radio.qoi <- logitsimfd(gabon.radio.scen, gabon.radio.simbetas, ci=0.95) #and get quantities of interest
gabon.radio.qoi.df <- as.data.frame(gabon.radio.qoi) %>% mutate(medium = 'radio', country = "gabon")

#Ghana
mdata <- extractdata(radio.country, ghana_data, na.rm = TRUE) #only extract data used in model
ghana.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
ghana.radio.pe <- fixef(ghana.radio.result) # extract pe
ghana.radio.vc <- vcov(ghana.radio.result) # extract vcov
ghana.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
ghana.radio.scen <- cfName(ghana.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
ghana.radio.scen <- cfChange(ghana.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
ghana.radio.simbetas <- mvrnorm(sims, ghana.radio.pe, ghana.radio.vc) #create simbetas
ghana.radio.qoi <- logitsimfd(ghana.radio.scen, ghana.radio.simbetas, ci=0.95) #and get quantities of interest
ghana.radio.qoi.df <- as.data.frame(ghana.radio.qoi) %>% mutate(medium = 'radio', country = "ghana")

#Guinea
mdata <- extractdata(radio.country, guinea_data, na.rm = TRUE) #only extract data used in model
guinea.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
guinea.radio.pe <- fixef(guinea.radio.result) # extract pe
guinea.radio.vc <- vcov(guinea.radio.result) # extract vcov
guinea.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
guinea.radio.scen <- cfName(guinea.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
guinea.radio.scen <- cfChange(guinea.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
guinea.radio.simbetas <- mvrnorm(sims, guinea.radio.pe, guinea.radio.vc) #create simbetas
guinea.radio.qoi <- logitsimfd(guinea.radio.scen, guinea.radio.simbetas, ci=0.95) #and get quantities of interest
guinea.radio.qoi.df <- as.data.frame(guinea.radio.qoi) %>% mutate(medium = 'radio', country = "guinea")

#Kenya
mdata <- extractdata(radio.country, kenya_data, na.rm = TRUE) #only extract data used in model
kenya.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
kenya.radio.pe <- fixef(kenya.radio.result) # extract pe
kenya.radio.vc <- vcov(kenya.radio.result) # extract vcov
kenya.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
kenya.radio.scen <- cfName(kenya.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
kenya.radio.scen <- cfChange(kenya.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
kenya.radio.simbetas <- mvrnorm(sims, kenya.radio.pe, kenya.radio.vc) #create simbetas
kenya.radio.qoi <- logitsimfd(kenya.radio.scen, kenya.radio.simbetas, ci=0.95) #and get quantities of interest
kenya.radio.qoi.df <- as.data.frame(kenya.radio.qoi) %>% mutate(medium = 'radio', country = "kenya")

#Lesotho
mdata <- extractdata(radio.country, lesotho_data, na.rm = TRUE) #only extract data used in model
lesotho.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
lesotho.radio.pe <- fixef(lesotho.radio.result) # extract pe
lesotho.radio.vc <- vcov(lesotho.radio.result) # extract vcov
lesotho.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
lesotho.radio.scen <- cfName(lesotho.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
lesotho.radio.scen <- cfChange(lesotho.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
lesotho.radio.simbetas <- mvrnorm(sims, lesotho.radio.pe, lesotho.radio.vc) #create simbetas
lesotho.radio.qoi <- logitsimfd(lesotho.radio.scen, lesotho.radio.simbetas, ci=0.95) #and get quantities of interest
lesotho.radio.qoi.df <- as.data.frame(lesotho.radio.qoi) %>% mutate(medium = 'radio', country = "lesotho")

#liberia
mdata <- extractdata(radio.country, liberia_data, na.rm = TRUE) #only extract data used in model
liberia.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
liberia.radio.pe <- fixef(liberia.radio.result) # extract pe
liberia.radio.vc <- vcov(liberia.radio.result) # extract vcov
liberia.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
liberia.radio.scen <- cfName(liberia.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
liberia.radio.scen <- cfChange(liberia.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
liberia.radio.simbetas <- mvrnorm(sims, liberia.radio.pe, liberia.radio.vc) #create simbetas
liberia.radio.qoi <- logitsimfd(liberia.radio.scen, liberia.radio.simbetas, ci=0.95) #and get quantities of interest
liberia.radio.qoi.df <- as.data.frame(liberia.radio.qoi) %>% mutate(medium = 'radio', country = "liberia")

#madagascar
mdata <- extractdata(radio.country, madag_data, na.rm = TRUE) #only extract data used in model
madag.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
madag.radio.pe <- fixef(madag.radio.result) # extract pe
madag.radio.vc <- vcov(madag.radio.result) # extract vcov
madag.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
madag.radio.scen <- cfName(madag.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
madag.radio.scen <- cfChange(madag.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
madag.radio.simbetas <- mvrnorm(sims, madag.radio.pe, madag.radio.vc) #create simbetas
madag.radio.qoi <- logitsimfd(madag.radio.scen, madag.radio.simbetas, ci=0.95) #and get quantities of interest
madag.radio.qoi.df <- as.data.frame(madag.radio.qoi) %>% mutate(medium = 'radio', country = "madagascar")

#malawi
mdata <- extractdata(radio.country, malawi_data, na.rm = TRUE) #only extract data used in model
malawi.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
malawi.radio.pe <- fixef(malawi.radio.result) # extract pe
malawi.radio.vc <- vcov(malawi.radio.result) # extract vcov
malawi.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
malawi.radio.scen <- cfName(malawi.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
malawi.radio.scen <- cfChange(malawi.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
malawi.radio.simbetas <- mvrnorm(sims, malawi.radio.pe, malawi.radio.vc) #create simbetas
malawi.radio.qoi <- logitsimfd(malawi.radio.scen, malawi.radio.simbetas, ci=0.95) #and get quantities of interest
malawi.radio.qoi.df <- as.data.frame(malawi.radio.qoi) %>% mutate(medium = 'radio', country = "malawi")

#mali
mdata <- extractdata(radio.country, mali_data, na.rm = TRUE) #only extract data used in model
mali.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
mali.radio.pe <- fixef(mali.radio.result) # extract pe
mali.radio.vc <- vcov(mali.radio.result) # extract vcov
mali.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mali.radio.scen <- cfName(mali.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
mali.radio.scen <- cfChange(mali.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
mali.radio.simbetas <- mvrnorm(sims, mali.radio.pe, mali.radio.vc) #create simbetas
mali.radio.qoi <- logitsimfd(mali.radio.scen, mali.radio.simbetas, ci=0.95) #and get quantities of interest
mali.radio.qoi.df <- as.data.frame(mali.radio.qoi) %>% mutate(medium = 'radio', country = "mali")

#mauritius
mdata <- extractdata(radio.country, maurit_data, na.rm = TRUE) #only extract data used in model
maurit.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
maurit.radio.pe <- fixef(maurit.radio.result) # extract pe
maurit.radio.vc <- vcov(maurit.radio.result) # extract vcov
maurit.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
maurit.radio.scen <- cfName(maurit.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
maurit.radio.scen <- cfChange(maurit.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
maurit.radio.simbetas <- mvrnorm(sims, maurit.radio.pe, maurit.radio.vc) #create simbetas
maurit.radio.qoi <- logitsimfd(maurit.radio.scen, maurit.radio.simbetas, ci=0.95) #and get quantities of interest
maurit.radio.qoi.df <- as.data.frame(maurit.radio.qoi) %>% mutate(medium = 'radio', country = "mauritius")

#morocco
mdata <- extractdata(radio.country, morocco_data, na.rm = TRUE) #only extract data used in model
morocco.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
morocco.radio.pe <- fixef(morocco.radio.result) # extract pe
morocco.radio.vc <- vcov(morocco.radio.result) # extract vcov
morocco.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
morocco.radio.scen <- cfName(morocco.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
morocco.radio.scen <- cfChange(morocco.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
morocco.radio.simbetas <- mvrnorm(sims, morocco.radio.pe, morocco.radio.vc) #create simbetas
morocco.radio.qoi <- logitsimfd(morocco.radio.scen, morocco.radio.simbetas, ci=0.95) #and get quantities of interest
morocco.radio.qoi.df <- as.data.frame(morocco.radio.qoi) %>% mutate(medium = 'radio', country = "morocco")

#mozambique
mdata <- extractdata(radio.country, mozamb_data, na.rm = TRUE) #only extract data used in model
mozamb.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
mozamb.radio.pe <- fixef(mozamb.radio.result) # extract pe
mozamb.radio.vc <- vcov(mozamb.radio.result) # extract vcov
mozamb.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mozamb.radio.scen <- cfName(mozamb.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
mozamb.radio.scen <- cfChange(mozamb.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
mozamb.radio.simbetas <- mvrnorm(sims, mozamb.radio.pe, mozamb.radio.vc) #create simbetas
mozamb.radio.qoi <- logitsimfd(mozamb.radio.scen, mozamb.radio.simbetas, ci=0.95) #and get quantities of interest
mozamb.radio.qoi.df <- as.data.frame(mozamb.radio.qoi) %>% mutate(medium = 'radio', country = "mozambique")

#namibia
mdata <- extractdata(radio.country, namibia_data, na.rm = TRUE) #only extract data used in model
namibia.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
namibia.radio.pe <- fixef(namibia.radio.result) # extract pe
namibia.radio.vc <- vcov(namibia.radio.result) # extract vcov
namibia.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
namibia.radio.scen <- cfName(namibia.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
namibia.radio.scen <- cfChange(namibia.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
namibia.radio.simbetas <- mvrnorm(sims, namibia.radio.pe, namibia.radio.vc) #create simbetas
namibia.radio.qoi <- logitsimfd(namibia.radio.scen, namibia.radio.simbetas, ci=0.95) #and get quantities of interest
namibia.radio.qoi.df <- as.data.frame(namibia.radio.qoi) %>% mutate(medium = 'radio', country = "namibia")

#niger
mdata <- extractdata(radio.country, niger_data, na.rm = TRUE) #only extract data used in model
niger.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
niger.radio.pe <- fixef(niger.radio.result) # extract pe
niger.radio.vc <- vcov(niger.radio.result) # extract vcov
niger.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
niger.radio.scen <- cfName(niger.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
niger.radio.scen <- cfChange(niger.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
niger.radio.simbetas <- mvrnorm(sims, niger.radio.pe, niger.radio.vc) #create simbetas
niger.radio.qoi <- logitsimfd(niger.radio.scen, niger.radio.simbetas, ci=0.95) #and get quantities of interest
niger.radio.qoi.df <- as.data.frame(niger.radio.qoi) %>% mutate(medium = 'radio', country = "niger")

#nigeria
mdata <- extractdata(radio.country, nigeria_data, na.rm = TRUE) #only extract data used in model
nigeria.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
nigeria.radio.pe <- fixef(nigeria.radio.result) # extract pe
nigeria.radio.vc <- vcov(nigeria.radio.result) # extract vcov
nigeria.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
nigeria.radio.scen <- cfName(nigeria.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
nigeria.radio.scen <- cfChange(nigeria.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
nigeria.radio.simbetas <- mvrnorm(sims, nigeria.radio.pe, nigeria.radio.vc) #create simbetas
nigeria.radio.qoi <- logitsimfd(nigeria.radio.scen, nigeria.radio.simbetas, ci=0.95) #and get quantities of interest
nigeria.radio.qoi.df <- as.data.frame(nigeria.radio.qoi) %>% mutate(medium = 'radio', country = "nigeria")

#Sao Tome Principe
mdata <- extractdata(radio.country, stp_data, na.rm = TRUE) #only extract data used in model
stp.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
stp.radio.pe <- fixef(stp.radio.result) # extract pe
stp.radio.vc <- vcov(stp.radio.result) # extract vcov
stp.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
stp.radio.scen <- cfName(stp.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
stp.radio.scen <- cfChange(stp.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
stp.radio.simbetas <- mvrnorm(sims, stp.radio.pe, stp.radio.vc) #create simbetas
stp.radio.qoi <- logitsimfd(stp.radio.scen, stp.radio.simbetas, ci=0.95) #and get quantities of interest
stp.radio.qoi.df <- as.data.frame(stp.radio.qoi) %>% mutate(medium = 'radio', country = "são tomé and príncipe")

#senegal
mdata <- extractdata(radio.country, senegal_data, na.rm = TRUE) #only extract data used in model
senegal.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
senegal.radio.pe <- fixef(senegal.radio.result) # extract pe
senegal.radio.vc <- vcov(senegal.radio.result) # extract vcov
senegal.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
senegal.radio.scen <- cfName(senegal.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
senegal.radio.scen <- cfChange(senegal.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
senegal.radio.simbetas <- mvrnorm(sims, senegal.radio.pe, senegal.radio.vc) #create simbetas
senegal.radio.qoi <- logitsimfd(senegal.radio.scen, senegal.radio.simbetas, ci=0.95) #and get quantities of interest
senegal.radio.qoi.df <- as.data.frame(senegal.radio.qoi) %>% mutate(medium = 'radio', country = "senegal")

#sierra leone
mdata <- extractdata(radio.country, sierra_data, na.rm = TRUE) #only extract data used in model
sierra.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
sierra.radio.pe <- fixef(sierra.radio.result) # extract pe
sierra.radio.vc <- vcov(sierra.radio.result) # extract vcov
sierra.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
sierra.radio.scen <- cfName(sierra.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
sierra.radio.scen <- cfChange(sierra.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
sierra.radio.simbetas <- mvrnorm(sims, sierra.radio.pe, sierra.radio.vc) #create simbetas
sierra.radio.qoi <- logitsimfd(sierra.radio.scen, sierra.radio.simbetas, ci=0.95) #and get quantities of interest
sierra.radio.qoi.df <- as.data.frame(sierra.radio.qoi) %>% mutate(medium = 'radio', country = "sierra leone")

#south africa
mdata <- extractdata(radio.country, southaf_data, na.rm = TRUE) #only extract data used in model
southaf.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
southaf.radio.pe <- fixef(southaf.radio.result) # extract pe
southaf.radio.vc <- vcov(southaf.radio.result) # extract vcov
southaf.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
southaf.radio.scen <- cfName(southaf.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
southaf.radio.scen <- cfChange(southaf.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
southaf.radio.simbetas <- mvrnorm(sims, southaf.radio.pe, southaf.radio.vc) #create simbetas
southaf.radio.qoi <- logitsimfd(southaf.radio.scen, southaf.radio.simbetas, ci=0.95) #and get quantities of interest
southaf.radio.qoi.df <- as.data.frame(southaf.radio.qoi) %>% mutate(medium = 'radio', country = "south africa")

#swaziland
mdata <- extractdata(radio.country, swaz_data, na.rm = TRUE) #only extract data used in model
swaz.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
swaz.radio.pe <- fixef(swaz.radio.result) # extract pe
swaz.radio.vc <- vcov(swaz.radio.result) # extract vcov
swaz.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
swaz.radio.scen <- cfName(swaz.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
swaz.radio.scen <- cfChange(swaz.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
swaz.radio.simbetas <- mvrnorm(sims, swaz.radio.pe, swaz.radio.vc) #create simbetas
swaz.radio.qoi <- logitsimfd(swaz.radio.scen, swaz.radio.simbetas, ci=0.95) #and get quantities of interest
swaz.radio.qoi.df <- as.data.frame(swaz.radio.qoi) %>% mutate(medium = 'radio', country = "swaziland")

#tanzania
mdata <- extractdata(radio.country, tanz_data, na.rm = TRUE) #only extract data used in model
tanz.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
tanz.radio.pe <- fixef(tanz.radio.result) # extract pe
tanz.radio.vc <- vcov(tanz.radio.result) # extract vcov
tanz.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tanz.radio.scen <- cfName(tanz.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
tanz.radio.scen <- cfChange(tanz.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
tanz.radio.simbetas <- mvrnorm(sims, tanz.radio.pe, tanz.radio.vc) #create simbetas
tanz.radio.qoi <- logitsimfd(tanz.radio.scen, tanz.radio.simbetas, ci=0.95) #and get quantities of interest
tanz.radio.qoi.df <- as.data.frame(tanz.radio.qoi) %>% mutate(medium = 'radio', country = "tanzania")

#togo
mdata <- extractdata(radio.country, togo_data, na.rm = TRUE) #only extract data used in model
togo.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
togo.radio.pe <- fixef(togo.radio.result) # extract pe
togo.radio.vc <- vcov(togo.radio.result) # extract vcov
togo.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
togo.radio.scen <- cfName(togo.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
togo.radio.scen <- cfChange(togo.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
togo.radio.simbetas <- mvrnorm(sims, togo.radio.pe, togo.radio.vc) #create simbetas
togo.radio.qoi <- logitsimfd(togo.radio.scen, togo.radio.simbetas, ci=0.95) #and get quantities of interest
togo.radio.qoi.df <- as.data.frame(togo.radio.qoi) %>% mutate(medium = 'radio', country = "togo")

#tunisia
mdata <- extractdata(radio.country, tunisia_data, na.rm = TRUE) #only extract data used in model
tunisia.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
tunisia.radio.pe <- fixef(tunisia.radio.result) # extract pe
tunisia.radio.vc <- vcov(tunisia.radio.result) # extract vcov
tunisia.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tunisia.radio.scen <- cfName(tunisia.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
tunisia.radio.scen <- cfChange(tunisia.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
tunisia.radio.simbetas <- mvrnorm(sims, tunisia.radio.pe, tunisia.radio.vc) #create simbetas
tunisia.radio.qoi <- logitsimfd(tunisia.radio.scen, tunisia.radio.simbetas, ci=0.95) #and get quantities of interest
tunisia.radio.qoi.df <- as.data.frame(tunisia.radio.qoi) %>% mutate(medium = 'radio', country = "tunisia")

#uganda
mdata <- extractdata(radio.country, uganda_data, na.rm = TRUE) #only extract data used in model
uganda.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
uganda.radio.pe <- fixef(uganda.radio.result) # extract pe
uganda.radio.vc <- vcov(uganda.radio.result) # extract vcov
uganda.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
uganda.radio.scen <- cfName(uganda.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
uganda.radio.scen <- cfChange(uganda.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
uganda.radio.simbetas <- mvrnorm(sims, uganda.radio.pe, uganda.radio.vc) #create simbetas
uganda.radio.qoi <- logitsimfd(uganda.radio.scen, uganda.radio.simbetas, ci=0.95) #and get quantities of interest
uganda.radio.qoi.df <- as.data.frame(uganda.radio.qoi) %>% mutate(medium = 'radio', country = "uganda")

#zambia
mdata <- extractdata(radio.country, zambia_data, na.rm = TRUE) #only extract data used in model
zambia.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
zambia.radio.pe <- fixef(zambia.radio.result) # extract pe
zambia.radio.vc <- vcov(zambia.radio.result) # extract vcov
zambia.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zambia.radio.scen <- cfName(zambia.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
zambia.radio.scen <- cfChange(zambia.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
zambia.radio.simbetas <- mvrnorm(sims, zambia.radio.pe, zambia.radio.vc) #create simbetas
zambia.radio.qoi <- logitsimfd(zambia.radio.scen, zambia.radio.simbetas, ci=0.95) #and get quantities of interest
zambia.radio.qoi.df <- as.data.frame(zambia.radio.qoi) %>% mutate(medium = 'radio', country = "zambia")

#zimbabwe
mdata <- extractdata(radio.country, zimbab_data, na.rm = TRUE) #only extract data used in model
zimbab.radio.result <- glmer(radio.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
zimbab.radio.pe <- fixef(zimbab.radio.result) # extract pe
zimbab.radio.vc <- vcov(zimbab.radio.result) # extract vcov
zimbab.radio.scen <- cfMake(radio.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zimbab.radio.scen <- cfName(zimbab.radio.scen, "radio 1 vs 5", scen=1) #configure scenario of 1 to 5
zimbab.radio.scen <- cfChange(zimbab.radio.scen, "radio", 
                               x = 5,
                               xpre=1,
                               scen=1)
zimbab.radio.simbetas <- mvrnorm(sims, zimbab.radio.pe, zimbab.radio.vc) #create simbetas
zimbab.radio.qoi <- logitsimfd(zimbab.radio.scen, zimbab.radio.simbetas, ci=0.95) #and get quantities of interest
zimbab.radio.qoi.df <- as.data.frame(zimbab.radio.qoi) %>% mutate(medium = 'radio', country = "zimbabwe")

##########
# TV
##########

tv.country <- #use this model to estimate result
  sexuality2 ~
  tv +
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 | district)  #add VI for District

tv.country.nodist <- #use this model to create scenarios (drop district)
  sexuality2 ~
  tv +
  media_notv + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban

#Benin
mdata <- extractdata(tv.country, benin_data, na.rm = TRUE) #only extract data used in model
benin.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                            nAGQ = 0 ) # optimizer for faster but less precise fit
benin.tv.pe <- fixef(benin.tv.result) # extract pe
benin.tv.vc <- vcov(benin.tv.result) # extract vcov
benin.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
benin.tv.scen <- cfName(benin.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
benin.tv.scen <- cfChange(benin.tv.scen, "tv", 
                             x = 5,
                             xpre=1,
                             scen=1)
benin.tv.simbetas <- mvrnorm(sims, benin.tv.pe, benin.tv.vc) #create simbetas
benin.tv.qoi <- logitsimfd(benin.tv.scen, benin.tv.simbetas, ci=0.95) #and get quantities of interest
benin.tv.qoi.df <- as.data.frame(benin.tv.qoi) %>% mutate(medium = 'tv', country = "benin")

#botswana
mdata <- extractdata(tv.country.nodist, botswana_data, na.rm = TRUE) #only extract data used in model
botswana.tv.result <- glm(tv.country.nodist, family=binomial, data=mdata) # estimate reg glm b/c no districts in botswana
botswana.tv.pe <- botswana.tv.result$coefficients # extract pe
botswana.tv.vc <- vcov(botswana.tv.result) # extract vcov
botswana.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
botswana.tv.scen <- cfName(botswana.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
botswana.tv.scen <- cfChange(botswana.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
botswana.tv.simbetas <- mvrnorm(sims, botswana.tv.pe, botswana.tv.vc) #create simbetas
botswana.tv.qoi <- logitsimfd(botswana.tv.scen, botswana.tv.simbetas, ci=0.95) #and get quantities of interest
botswana.tv.qoi.df <- as.data.frame(botswana.tv.qoi) %>% mutate(medium = 'tv', country = "botswana")

#burkina faso
mdata <- extractdata(tv.country, burkina_data, na.rm = TRUE) #only extract data used in model
burkina.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
burkina.tv.pe <- fixef(burkina.tv.result) # extract pe
burkina.tv.vc <- vcov(burkina.tv.result) # extract vcov
burkina.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burkina.tv.scen <- cfName(burkina.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
burkina.tv.scen <- cfChange(burkina.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
burkina.tv.simbetas <- mvrnorm(sims, burkina.tv.pe, burkina.tv.vc) #create simbetas
burkina.tv.qoi <- logitsimfd(burkina.tv.scen, burkina.tv.simbetas, ci=0.95) #and get quantities of interest
burkina.tv.qoi.df <- as.data.frame(burkina.tv.qoi) %>% mutate(medium = 'tv', country = "burkina faso")

#burundi
mdata <- extractdata(tv.country, burundi_data, na.rm = TRUE) #only extract data used in model
burundi.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
burundi.tv.pe <- fixef(burundi.tv.result) # extract pe
burundi.tv.vc <- vcov(burundi.tv.result) # extract vcov
burundi.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burundi.tv.scen <- cfName(burundi.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
burundi.tv.scen <- cfChange(burundi.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
burundi.tv.simbetas <- mvrnorm(sims, burundi.tv.pe, burundi.tv.vc) #create simbetas
burundi.tv.qoi <- logitsimfd(burundi.tv.scen, burundi.tv.simbetas, ci=0.95) #and get quantities of interest
burundi.tv.qoi.df <- as.data.frame(burundi.tv.qoi) %>% mutate(medium = 'tv', country = "burundi")

#cameroon
mdata <- extractdata(tv.country, cameroon_data, na.rm = TRUE) #only extract data used in model
cameroon.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
cameroon.tv.pe <- fixef(cameroon.tv.result) # extract pe
cameroon.tv.vc <- vcov(cameroon.tv.result) # extract vcov
cameroon.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cameroon.tv.scen <- cfName(cameroon.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
cameroon.tv.scen <- cfChange(cameroon.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
cameroon.tv.simbetas <- mvrnorm(sims, cameroon.tv.pe, cameroon.tv.vc) #create simbetas
cameroon.tv.qoi <- logitsimfd(cameroon.tv.scen, cameroon.tv.simbetas, ci=0.95) #and get quantities of interest
cameroon.tv.qoi.df <- as.data.frame(cameroon.tv.qoi) %>% mutate(medium = 'tv', country = "cameroon")

#cape verde
mdata <- extractdata(tv.country, capeverde_data, na.rm = TRUE) #only extract data used in model
capeverde.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
capeverde.tv.pe <- fixef(capeverde.tv.result) # extract pe
capeverde.tv.vc <- vcov(capeverde.tv.result) # extract vcov
capeverde.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
capeverde.tv.scen <- cfName(capeverde.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
capeverde.tv.scen <- cfChange(capeverde.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
capeverde.tv.simbetas <- mvrnorm(sims, capeverde.tv.pe, capeverde.tv.vc) #create simbetas
capeverde.tv.qoi <- logitsimfd(capeverde.tv.scen, capeverde.tv.simbetas, ci=0.95) #and get quantities of interest
capeverde.tv.qoi.df <- as.data.frame(capeverde.tv.qoi) %>% mutate(medium = 'tv', country = "cape verde")

#cote d'ivoire
mdata <- extractdata(tv.country, cdi_data, na.rm = TRUE) #only extract data used in model
cdi.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
cdi.tv.pe <- fixef(cdi.tv.result) # extract pe
cdi.tv.vc <- vcov(cdi.tv.result) # extract vcov
cdi.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cdi.tv.scen <- cfName(cdi.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
cdi.tv.scen <- cfChange(cdi.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
cdi.tv.simbetas <- mvrnorm(sims, cdi.tv.pe, cdi.tv.vc) #create simbetas
cdi.tv.qoi <- logitsimfd(cdi.tv.scen, cdi.tv.simbetas, ci=0.95) #and get quantities of interest
cdi.tv.qoi.df <- as.data.frame(cdi.tv.qoi) %>% mutate(medium = 'tv', country = "cote d'ivoire")

#gabon
mdata <- extractdata(tv.country, gabon_data, na.rm = TRUE) #only extract data used in model
gabon.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
gabon.tv.pe <- fixef(gabon.tv.result) # extract pe
gabon.tv.vc <- vcov(gabon.tv.result) # extract vcov
gabon.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
gabon.tv.scen <- cfName(gabon.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
gabon.tv.scen <- cfChange(gabon.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
gabon.tv.simbetas <- mvrnorm(sims, gabon.tv.pe, gabon.tv.vc) #create simbetas
gabon.tv.qoi <- logitsimfd(gabon.tv.scen, gabon.tv.simbetas, ci=0.95) #and get quantities of interest
gabon.tv.qoi.df <- as.data.frame(gabon.tv.qoi) %>% mutate(medium = 'tv', country = "gabon")

#ghana
mdata <- extractdata(tv.country, ghana_data, na.rm = TRUE) #only extract data used in model
ghana.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
ghana.tv.pe <- fixef(ghana.tv.result) # extract pe
ghana.tv.vc <- vcov(ghana.tv.result) # extract vcov
ghana.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
ghana.tv.scen <- cfName(ghana.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
ghana.tv.scen <- cfChange(ghana.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
ghana.tv.simbetas <- mvrnorm(sims, ghana.tv.pe, ghana.tv.vc) #create simbetas
ghana.tv.qoi <- logitsimfd(ghana.tv.scen, ghana.tv.simbetas, ci=0.95) #and get quantities of interest
ghana.tv.qoi.df <- as.data.frame(ghana.tv.qoi) %>% mutate(medium = 'tv', country = "ghana")

#guinea 
mdata <- extractdata(tv.country, guinea_data, na.rm = TRUE) #only extract data used in model
guinea.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
guinea.tv.pe <- fixef(guinea.tv.result) # extract pe
guinea.tv.vc <- vcov(guinea.tv.result) # extract vcov
guinea.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
guinea.tv.scen <- cfName(guinea.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
guinea.tv.scen <- cfChange(guinea.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
guinea.tv.simbetas <- mvrnorm(sims, guinea.tv.pe, guinea.tv.vc) #create simbetas
guinea.tv.qoi <- logitsimfd(guinea.tv.scen, guinea.tv.simbetas, ci=0.95) #and get quantities of interest
guinea.tv.qoi.df <- as.data.frame(guinea.tv.qoi) %>% mutate(medium = 'tv', country = "guinea")

#kenya
mdata <- extractdata(tv.country, kenya_data, na.rm = TRUE) #only extract data used in model
kenya.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
kenya.tv.pe <- fixef(kenya.tv.result) # extract pe
kenya.tv.vc <- vcov(kenya.tv.result) # extract vcov
kenya.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
kenya.tv.scen <- cfName(kenya.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
kenya.tv.scen <- cfChange(kenya.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
kenya.tv.simbetas <- mvrnorm(sims, kenya.tv.pe, kenya.tv.vc) #create simbetas
kenya.tv.qoi <- logitsimfd(kenya.tv.scen, kenya.tv.simbetas, ci=0.95) #and get quantities of interest
kenya.tv.qoi.df <- as.data.frame(kenya.tv.qoi) %>% mutate(medium = 'tv', country = "kenya")

#lesotho
mdata <- extractdata(tv.country, lesotho_data, na.rm = TRUE) #only extract data used in model
lesotho.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
lesotho.tv.pe <- fixef(lesotho.tv.result) # extract pe
lesotho.tv.vc <- vcov(lesotho.tv.result) # extract vcov
lesotho.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
lesotho.tv.scen <- cfName(lesotho.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
lesotho.tv.scen <- cfChange(lesotho.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
lesotho.tv.simbetas <- mvrnorm(sims, lesotho.tv.pe, lesotho.tv.vc) #create simbetas
lesotho.tv.qoi <- logitsimfd(lesotho.tv.scen, lesotho.tv.simbetas, ci=0.95) #and get quantities of interest
lesotho.tv.qoi.df <- as.data.frame(lesotho.tv.qoi) %>% mutate(medium = 'tv', country = "lesotho")

#liberia
mdata <- extractdata(tv.country, liberia_data, na.rm = TRUE) #only extract data used in model
liberia.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
liberia.tv.pe <- fixef(liberia.tv.result) # extract pe
liberia.tv.vc <- vcov(liberia.tv.result) # extract vcov
liberia.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
liberia.tv.scen <- cfName(liberia.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
liberia.tv.scen <- cfChange(liberia.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
liberia.tv.simbetas <- mvrnorm(sims, liberia.tv.pe, liberia.tv.vc) #create simbetas
liberia.tv.qoi <- logitsimfd(liberia.tv.scen, liberia.tv.simbetas, ci=0.95) #and get quantities of interest
liberia.tv.qoi.df <- as.data.frame(liberia.tv.qoi) %>% mutate(medium = 'tv', country = "liberia")

#madagascar
mdata <- extractdata(tv.country, madag_data, na.rm = TRUE) #only extract data used in model
madag.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
madag.tv.pe <- fixef(madag.tv.result) # extract pe
madag.tv.vc <- vcov(madag.tv.result) # extract vcov
madag.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
madag.tv.scen <- cfName(madag.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
madag.tv.scen <- cfChange(madag.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
madag.tv.simbetas <- mvrnorm(sims, madag.tv.pe, madag.tv.vc) #create simbetas
madag.tv.qoi <- logitsimfd(madag.tv.scen, madag.tv.simbetas, ci=0.95) #and get quantities of interest
madag.tv.qoi.df <- as.data.frame(madag.tv.qoi) %>% mutate(medium = 'tv', country = "madagascar")

#malawi
mdata <- extractdata(tv.country, malawi_data, na.rm = TRUE) #only extract data used in model
malawi.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
malawi.tv.pe <- fixef(malawi.tv.result) # extract pe
malawi.tv.vc <- vcov(malawi.tv.result) # extract vcov
malawi.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
malawi.tv.scen <- cfName(malawi.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
malawi.tv.scen <- cfChange(malawi.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
malawi.tv.simbetas <- mvrnorm(sims, malawi.tv.pe, malawi.tv.vc) #create simbetas
malawi.tv.qoi <- logitsimfd(malawi.tv.scen, malawi.tv.simbetas, ci=0.95) #and get quantities of interest
malawi.tv.qoi.df <- as.data.frame(malawi.tv.qoi) %>% mutate(medium = 'tv', country = "malawi")

#mali
mdata <- extractdata(tv.country, mali_data, na.rm = TRUE) #only extract data used in model
mali.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
mali.tv.pe <- fixef(mali.tv.result) # extract pe
mali.tv.vc <- vcov(mali.tv.result) # extract vcov
mali.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mali.tv.scen <- cfName(mali.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
mali.tv.scen <- cfChange(mali.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
mali.tv.simbetas <- mvrnorm(sims, mali.tv.pe, mali.tv.vc) #create simbetas
mali.tv.qoi <- logitsimfd(mali.tv.scen, mali.tv.simbetas, ci=0.95) #and get quantities of interest
mali.tv.qoi.df <- as.data.frame(mali.tv.qoi) %>% mutate(medium = 'tv', country = "mali")

#mauritius
mdata <- extractdata(tv.country, maurit_data, na.rm = TRUE) #only extract data used in model
maurit.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
maurit.tv.pe <- fixef(maurit.tv.result) # extract pe
maurit.tv.vc <- vcov(maurit.tv.result) # extract vcov
maurit.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
maurit.tv.scen <- cfName(maurit.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
maurit.tv.scen <- cfChange(maurit.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
maurit.tv.simbetas <- mvrnorm(sims, maurit.tv.pe, maurit.tv.vc) #create simbetas
maurit.tv.qoi <- logitsimfd(maurit.tv.scen, maurit.tv.simbetas, ci=0.95) #and get quantities of interest
maurit.tv.qoi.df <- as.data.frame(maurit.tv.qoi) %>% mutate(medium = 'tv', country = "mauritius")

#morocco
mdata <- extractdata(tv.country, morocco_data, na.rm = TRUE) #only extract data used in model
morocco.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
morocco.tv.pe <- fixef(morocco.tv.result) # extract pe
morocco.tv.vc <- vcov(morocco.tv.result) # extract vcov
morocco.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
morocco.tv.scen <- cfName(morocco.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
morocco.tv.scen <- cfChange(morocco.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
morocco.tv.simbetas <- mvrnorm(sims, morocco.tv.pe, morocco.tv.vc) #create simbetas
morocco.tv.qoi <- logitsimfd(morocco.tv.scen, morocco.tv.simbetas, ci=0.95) #and get quantities of interest
morocco.tv.qoi.df <- as.data.frame(morocco.tv.qoi) %>% mutate(medium = 'tv', country = "morocco")

#mozambique
mdata <- extractdata(tv.country, mozamb_data, na.rm = TRUE) #only extract data used in model
mozamb.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
mozamb.tv.pe <- fixef(mozamb.tv.result) # extract pe
mozamb.tv.vc <- vcov(mozamb.tv.result) # extract vcov
mozamb.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mozamb.tv.scen <- cfName(mozamb.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
mozamb.tv.scen <- cfChange(mozamb.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
mozamb.tv.simbetas <- mvrnorm(sims, mozamb.tv.pe, mozamb.tv.vc) #create simbetas
mozamb.tv.qoi <- logitsimfd(mozamb.tv.scen, mozamb.tv.simbetas, ci=0.95) #and get quantities of interest
mozamb.tv.qoi.df <- as.data.frame(mozamb.tv.qoi) %>% mutate(medium = 'tv', country = "mozambique")

#namibia
mdata <- extractdata(tv.country, namibia_data, na.rm = TRUE) #only extract data used in model
namibia.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
namibia.tv.pe <- fixef(namibia.tv.result) # extract pe
namibia.tv.vc <- vcov(namibia.tv.result) # extract vcov
namibia.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
namibia.tv.scen <- cfName(namibia.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
namibia.tv.scen <- cfChange(namibia.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
namibia.tv.simbetas <- mvrnorm(sims, namibia.tv.pe, namibia.tv.vc) #create simbetas
namibia.tv.qoi <- logitsimfd(namibia.tv.scen, namibia.tv.simbetas, ci=0.95) #and get quantities of interest
namibia.tv.qoi.df <- as.data.frame(namibia.tv.qoi) %>% mutate(medium = 'tv', country = "namibia")

#niger
mdata <- extractdata(tv.country, niger_data, na.rm = TRUE) #only extract data used in model
niger.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
niger.tv.pe <- fixef(niger.tv.result) # extract pe
niger.tv.vc <- vcov(niger.tv.result) # extract vcov
niger.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
niger.tv.scen <- cfName(niger.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
niger.tv.scen <- cfChange(niger.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
niger.tv.simbetas <- mvrnorm(sims, niger.tv.pe, niger.tv.vc) #create simbetas
niger.tv.qoi <- logitsimfd(niger.tv.scen, niger.tv.simbetas, ci=0.95) #and get quantities of interest
niger.tv.qoi.df <- as.data.frame(niger.tv.qoi) %>% mutate(medium = 'tv', country = "niger")

#nigeria
mdata <- extractdata(tv.country, nigeria_data, na.rm = TRUE) #only extract data used in model
nigeria.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
nigeria.tv.pe <- fixef(nigeria.tv.result) # extract pe
nigeria.tv.vc <- vcov(nigeria.tv.result) # extract vcov
nigeria.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
nigeria.tv.scen <- cfName(nigeria.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
nigeria.tv.scen <- cfChange(nigeria.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
nigeria.tv.simbetas <- mvrnorm(sims, nigeria.tv.pe, nigeria.tv.vc) #create simbetas
nigeria.tv.qoi <- logitsimfd(nigeria.tv.scen, nigeria.tv.simbetas, ci=0.95) #and get quantities of interest
nigeria.tv.qoi.df <- as.data.frame(nigeria.tv.qoi) %>% mutate(medium = 'tv', country = "nigeria")

#são tomé and príncipe
mdata <- extractdata(tv.country, stp_data, na.rm = TRUE) #only extract data used in model
stp.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
stp.tv.pe <- fixef(stp.tv.result) # extract pe
stp.tv.vc <- vcov(stp.tv.result) # extract vcov
stp.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
stp.tv.scen <- cfName(stp.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
stp.tv.scen <- cfChange(stp.tv.scen, "tv", 
                          x = 5,
                          xpre=1,
                          scen=1)
stp.tv.simbetas <- mvrnorm(sims, stp.tv.pe, stp.tv.vc) #create simbetas
stp.tv.qoi <- logitsimfd(stp.tv.scen, stp.tv.simbetas, ci=0.95) #and get quantities of interest
stp.tv.qoi.df <- as.data.frame(stp.tv.qoi) %>% mutate(medium = 'tv', country = "são tomé and príncipe")

#senegal
mdata <- extractdata(tv.country, senegal_data, na.rm = TRUE) #only extract data used in model
senegal.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
senegal.tv.pe <- fixef(senegal.tv.result) # extract pe
senegal.tv.vc <- vcov(senegal.tv.result) # extract vcov
senegal.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
senegal.tv.scen <- cfName(senegal.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
senegal.tv.scen <- cfChange(senegal.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
senegal.tv.simbetas <- mvrnorm(sims, senegal.tv.pe, senegal.tv.vc) #create simbetas
senegal.tv.qoi <- logitsimfd(senegal.tv.scen, senegal.tv.simbetas, ci=0.95) #and get quantities of interest
senegal.tv.qoi.df <- as.data.frame(senegal.tv.qoi) %>% mutate(medium = 'tv', country = "senegal")

#sierra
mdata <- extractdata(tv.country, sierra_data, na.rm = TRUE) #only extract data used in model
sierra.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
sierra.tv.pe <- fixef(sierra.tv.result) # extract pe
sierra.tv.vc <- vcov(sierra.tv.result) # extract vcov
sierra.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
sierra.tv.scen <- cfName(sierra.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
sierra.tv.scen <- cfChange(sierra.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
sierra.tv.simbetas <- mvrnorm(sims, sierra.tv.pe, sierra.tv.vc) #create simbetas
sierra.tv.qoi <- logitsimfd(sierra.tv.scen, sierra.tv.simbetas, ci=0.95) #and get quantities of interest
sierra.tv.qoi.df <- as.data.frame(sierra.tv.qoi) %>% mutate(medium = 'tv', country = "sierra leone")

#southaf
mdata <- extractdata(tv.country, southaf_data, na.rm = TRUE) #only extract data used in model
southaf.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
southaf.tv.pe <- fixef(southaf.tv.result) # extract pe
southaf.tv.vc <- vcov(southaf.tv.result) # extract vcov
southaf.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
southaf.tv.scen <- cfName(southaf.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
southaf.tv.scen <- cfChange(southaf.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
southaf.tv.simbetas <- mvrnorm(sims, southaf.tv.pe, southaf.tv.vc) #create simbetas
southaf.tv.qoi <- logitsimfd(southaf.tv.scen, southaf.tv.simbetas, ci=0.95) #and get quantities of interest
southaf.tv.qoi.df <- as.data.frame(southaf.tv.qoi) %>% mutate(medium = 'tv', country = "south africa")

#swaz
mdata <- extractdata(tv.country, swaz_data, na.rm = TRUE) #only extract data used in model
swaz.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
swaz.tv.pe <- fixef(swaz.tv.result) # extract pe
swaz.tv.vc <- vcov(swaz.tv.result) # extract vcov
swaz.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
swaz.tv.scen <- cfName(swaz.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
swaz.tv.scen <- cfChange(swaz.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
swaz.tv.simbetas <- mvrnorm(sims, swaz.tv.pe, swaz.tv.vc) #create simbetas
swaz.tv.qoi <- logitsimfd(swaz.tv.scen, swaz.tv.simbetas, ci=0.95) #and get quantities of interest
swaz.tv.qoi.df <- as.data.frame(swaz.tv.qoi) %>% mutate(medium = 'tv', country = "swaziland")

#tanz
mdata <- extractdata(tv.country, tanz_data, na.rm = TRUE) #only extract data used in model
tanz.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
tanz.tv.pe <- fixef(tanz.tv.result) # extract pe
tanz.tv.vc <- vcov(tanz.tv.result) # extract vcov
tanz.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tanz.tv.scen <- cfName(tanz.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
tanz.tv.scen <- cfChange(tanz.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
tanz.tv.simbetas <- mvrnorm(sims, tanz.tv.pe, tanz.tv.vc) #create simbetas
tanz.tv.qoi <- logitsimfd(tanz.tv.scen, tanz.tv.simbetas, ci=0.95) #and get quantities of interest
tanz.tv.qoi.df <- as.data.frame(tanz.tv.qoi) %>% mutate(medium = 'tv', country = "tanzania")

#togo
mdata <- extractdata(tv.country, togo_data, na.rm = TRUE) #only extract data used in model
togo.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
togo.tv.pe <- fixef(togo.tv.result) # extract pe
togo.tv.vc <- vcov(togo.tv.result) # extract vcov
togo.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
togo.tv.scen <- cfName(togo.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
togo.tv.scen <- cfChange(togo.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
togo.tv.simbetas <- mvrnorm(sims, togo.tv.pe, togo.tv.vc) #create simbetas
togo.tv.qoi <- logitsimfd(togo.tv.scen, togo.tv.simbetas, ci=0.95) #and get quantities of interest
togo.tv.qoi.df <- as.data.frame(togo.tv.qoi) %>% mutate(medium = 'tv', country = "togo")

#tunisia
mdata <- extractdata(tv.country, tunisia_data, na.rm = TRUE) #only extract data used in model
tunisia.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
tunisia.tv.pe <- fixef(tunisia.tv.result) # extract pe
tunisia.tv.vc <- vcov(tunisia.tv.result) # extract vcov
tunisia.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tunisia.tv.scen <- cfName(tunisia.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
tunisia.tv.scen <- cfChange(tunisia.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
tunisia.tv.simbetas <- mvrnorm(sims, tunisia.tv.pe, tunisia.tv.vc) #create simbetas
tunisia.tv.qoi <- logitsimfd(tunisia.tv.scen, tunisia.tv.simbetas, ci=0.95) #and get quantities of interest
tunisia.tv.qoi.df <- as.data.frame(tunisia.tv.qoi) %>% mutate(medium = 'tv', country = "tunisia")

#uganda
mdata <- extractdata(tv.country, uganda_data, na.rm = TRUE) #only extract data used in model
uganda.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
uganda.tv.pe <- fixef(uganda.tv.result) # extract pe
uganda.tv.vc <- vcov(uganda.tv.result) # extract vcov
uganda.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
uganda.tv.scen <- cfName(uganda.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
uganda.tv.scen <- cfChange(uganda.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
uganda.tv.simbetas <- mvrnorm(sims, uganda.tv.pe, uganda.tv.vc) #create simbetas
uganda.tv.qoi <- logitsimfd(uganda.tv.scen, uganda.tv.simbetas, ci=0.95) #and get quantities of interest
uganda.tv.qoi.df <- as.data.frame(uganda.tv.qoi) %>% mutate(medium = 'tv', country = "uganda")

#zambia
mdata <- extractdata(tv.country, zambia_data, na.rm = TRUE) #only extract data used in model
zambia.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
zambia.tv.pe <- fixef(zambia.tv.result) # extract pe
zambia.tv.vc <- vcov(zambia.tv.result) # extract vcov
zambia.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zambia.tv.scen <- cfName(zambia.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
zambia.tv.scen <- cfChange(zambia.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
zambia.tv.simbetas <- mvrnorm(sims, zambia.tv.pe, zambia.tv.vc) #create simbetas
zambia.tv.qoi <- logitsimfd(zambia.tv.scen, zambia.tv.simbetas, ci=0.95) #and get quantities of interest
zambia.tv.qoi.df <- as.data.frame(zambia.tv.qoi) %>% mutate(medium = 'tv', country = "zambia")

#zimbab
mdata <- extractdata(tv.country, zimbab_data, na.rm = TRUE) #only extract data used in model
zimbab.tv.result <- glmer(tv.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
zimbab.tv.pe <- fixef(zimbab.tv.result) # extract pe
zimbab.tv.vc <- vcov(zimbab.tv.result) # extract vcov
zimbab.tv.scen <- cfMake(tv.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zimbab.tv.scen <- cfName(zimbab.tv.scen, "tv 1 vs 5", scen=1) #configure scenario of 1 to 5
zimbab.tv.scen <- cfChange(zimbab.tv.scen, "tv", 
                            x = 5,
                            xpre=1,
                            scen=1)
zimbab.tv.simbetas <- mvrnorm(sims, zimbab.tv.pe, zimbab.tv.vc) #create simbetas
zimbab.tv.qoi <- logitsimfd(zimbab.tv.scen, zimbab.tv.simbetas, ci=0.95) #and get quantities of interest
zimbab.tv.qoi.df <- as.data.frame(zimbab.tv.qoi) %>% mutate(medium = 'tv', country = "zimbabwe")

############
# Newspaper
############

newspaper.country <- #use this model to estimate result
  sexuality2 ~
  newspaper +
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 | district)  #add VI for District

newspaper.country.nodist <- #use this model to create scenarios (drop district)
  sexuality2 ~
  newspaper +
  media_nopaper + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban

#Benin
mdata <- extractdata(newspaper.country, benin_data, na.rm = TRUE) #only extract data used in model
benin.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
benin.newspaper.pe <- fixef(benin.newspaper.result) # extract pe
benin.newspaper.vc <- vcov(benin.newspaper.result) # extract vcov
benin.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
benin.newspaper.scen <- cfName(benin.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
benin.newspaper.scen <- cfChange(benin.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
benin.newspaper.simbetas <- mvrnorm(sims, benin.newspaper.pe, benin.newspaper.vc) #create simbetas
benin.newspaper.qoi <- logitsimfd(benin.newspaper.scen, benin.newspaper.simbetas, ci=0.95) #and get quantities of interest
benin.newspaper.qoi.df <- as.data.frame(benin.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "benin")

#botswana
mdata <- extractdata(newspaper.country.nodist, botswana_data, na.rm = TRUE) #only extract data used in model
botswana.newspaper.result <- glm(newspaper.country.nodist, family=binomial, data=mdata)#estimate the model
botswana.newspaper.pe <- botswana.newspaper.result$coefficients # extract pe
botswana.newspaper.vc <- vcov(botswana.newspaper.result) # extract vcov
botswana.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
botswana.newspaper.scen <- cfName(botswana.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
botswana.newspaper.scen <- cfChange(botswana.newspaper.scen, "newspaper", 
                                x = 5,
                                xpre=1,
                                scen=1)
botswana.newspaper.simbetas <- mvrnorm(sims, botswana.newspaper.pe, botswana.newspaper.vc) #create simbetas
botswana.newspaper.qoi <- logitsimfd(botswana.newspaper.scen, botswana.newspaper.simbetas, ci=0.95) #and get quantities of interest
botswana.newspaper.qoi.df <- as.data.frame(botswana.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "botswana") #create df from qoi

#burkina faso
mdata <- extractdata(newspaper.country, burkina_data, na.rm = TRUE) #only extract data used in model
burkina.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
burkina.newspaper.pe <- fixef(burkina.newspaper.result) # extract pe
burkina.newspaper.vc <- vcov(burkina.newspaper.result) # extract vcov
burkina.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burkina.newspaper.scen <- cfName(burkina.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
burkina.newspaper.scen <- cfChange(burkina.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
burkina.newspaper.simbetas <- mvrnorm(sims, burkina.newspaper.pe, burkina.newspaper.vc) #create simbetas
burkina.newspaper.qoi <- logitsimfd(burkina.newspaper.scen, burkina.newspaper.simbetas, ci=0.95) #and get quantities of interest
burkina.newspaper.qoi.df <- as.data.frame(burkina.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "burkina faso")

#burundi
mdata <- extractdata(newspaper.country, burundi_data, na.rm = TRUE) #only extract data used in model
burundi.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
burundi.newspaper.pe <- fixef(burundi.newspaper.result) # extract pe
burundi.newspaper.vc <- vcov(burundi.newspaper.result) # extract vcov
burundi.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burundi.newspaper.scen <- cfName(burundi.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
burundi.newspaper.scen <- cfChange(burundi.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
burundi.newspaper.simbetas <- mvrnorm(sims, burundi.newspaper.pe, burundi.newspaper.vc) #create simbetas
burundi.newspaper.qoi <- logitsimfd(burundi.newspaper.scen, burundi.newspaper.simbetas, ci=0.95) #and get quantities of interest
burundi.newspaper.qoi.df <- as.data.frame(burundi.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "burundi")

#cameroon
mdata <- extractdata(newspaper.country, cameroon_data, na.rm = TRUE) #only extract data used in model
cameroon.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                            nAGQ = 0 ) # optimizer for faster but less precise fit
cameroon.newspaper.pe <- fixef(cameroon.newspaper.result) # extract pe
cameroon.newspaper.vc <- vcov(cameroon.newspaper.result) # extract vcov
cameroon.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cameroon.newspaper.scen <- cfName(cameroon.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
cameroon.newspaper.scen <- cfChange(cameroon.newspaper.scen, "newspaper", 
                             x = 5,
                             xpre=1,
                             scen=1)
cameroon.newspaper.simbetas <- mvrnorm(sims, cameroon.newspaper.pe, cameroon.newspaper.vc) #create simbetas
cameroon.newspaper.qoi <- logitsimfd(cameroon.newspaper.scen, cameroon.newspaper.simbetas, ci=0.95) #and get quantities of interest
cameroon.newspaper.qoi.df <- as.data.frame(cameroon.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "cameroon")

#cape verde
mdata <- extractdata(newspaper.country, capeverde_data, na.rm = TRUE) #only extract data used in model
capeverde.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                             nAGQ = 0 ) # optimizer for faster but less precise fit
capeverde.newspaper.pe <- fixef(capeverde.newspaper.result) # extract pe
capeverde.newspaper.vc <- vcov(capeverde.newspaper.result) # extract vcov
capeverde.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
capeverde.newspaper.scen <- cfName(capeverde.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
capeverde.newspaper.scen <- cfChange(capeverde.newspaper.scen, "newspaper", 
                              x = 5,
                              xpre=1,
                              scen=1)
capeverde.newspaper.simbetas <- mvrnorm(sims, capeverde.newspaper.pe, capeverde.newspaper.vc) #create simbetas
capeverde.newspaper.qoi <- logitsimfd(capeverde.newspaper.scen, capeverde.newspaper.simbetas, ci=0.95) #and get quantities of interest
capeverde.newspaper.qoi.df <- as.data.frame(capeverde.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "cape verde")

#cote d'ivoire
mdata <- extractdata(newspaper.country, cdi_data, na.rm = TRUE) #only extract data used in model
cdi.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                       nAGQ = 0 ) # optimizer for faster but less precise fit
cdi.newspaper.pe <- fixef(cdi.newspaper.result) # extract pe
cdi.newspaper.vc <- vcov(cdi.newspaper.result) # extract vcov
cdi.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cdi.newspaper.scen <- cfName(cdi.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
cdi.newspaper.scen <- cfChange(cdi.newspaper.scen, "newspaper", 
                        x = 5,
                        xpre=1,
                        scen=1)
cdi.newspaper.simbetas <- mvrnorm(sims, cdi.newspaper.pe, cdi.newspaper.vc) #create simbetas
cdi.newspaper.qoi <- logitsimfd(cdi.newspaper.scen, cdi.newspaper.simbetas, ci=0.95) #and get quantities of interest
cdi.newspaper.qoi.df <- as.data.frame(cdi.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "cote d'ivoire")

#gabon
mdata <- extractdata(newspaper.country, gabon_data, na.rm = TRUE) #only extract data used in model
gabon.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
gabon.newspaper.pe <- fixef(gabon.newspaper.result) # extract pe
gabon.newspaper.vc <- vcov(gabon.newspaper.result) # extract vcov
gabon.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
gabon.newspaper.scen <- cfName(gabon.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
gabon.newspaper.scen <- cfChange(gabon.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
gabon.newspaper.simbetas <- mvrnorm(sims, gabon.newspaper.pe, gabon.newspaper.vc) #create simbetas
gabon.newspaper.qoi <- logitsimfd(gabon.newspaper.scen, gabon.newspaper.simbetas, ci=0.95) #and get quantities of interest
gabon.newspaper.qoi.df <- as.data.frame(gabon.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "gabon")

#ghana
mdata <- extractdata(newspaper.country, ghana_data, na.rm = TRUE) #only extract data used in model
ghana.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
ghana.newspaper.pe <- fixef(ghana.newspaper.result) # extract pe
ghana.newspaper.vc <- vcov(ghana.newspaper.result) # extract vcov
ghana.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
ghana.newspaper.scen <- cfName(ghana.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
ghana.newspaper.scen <- cfChange(ghana.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
ghana.newspaper.simbetas <- mvrnorm(sims, ghana.newspaper.pe, ghana.newspaper.vc) #create simbetas
ghana.newspaper.qoi <- logitsimfd(ghana.newspaper.scen, ghana.newspaper.simbetas, ci=0.95) #and get quantities of interest
ghana.newspaper.qoi.df <- as.data.frame(ghana.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "ghana")

#guinea 
mdata <- extractdata(newspaper.country, guinea_data, na.rm = TRUE) #only extract data used in model
guinea.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
guinea.newspaper.pe <- fixef(guinea.newspaper.result) # extract pe
guinea.newspaper.vc <- vcov(guinea.newspaper.result) # extract vcov
guinea.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
guinea.newspaper.scen <- cfName(guinea.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
guinea.newspaper.scen <- cfChange(guinea.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
guinea.newspaper.simbetas <- mvrnorm(sims, guinea.newspaper.pe, guinea.newspaper.vc) #create simbetas
guinea.newspaper.qoi <- logitsimfd(guinea.newspaper.scen, guinea.newspaper.simbetas, ci=0.95) #and get quantities of interest
guinea.newspaper.qoi.df <- as.data.frame(guinea.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "guinea")

#kenya
mdata <- extractdata(newspaper.country, kenya_data, na.rm = TRUE) #only extract data used in model
kenya.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
kenya.newspaper.pe <- fixef(kenya.newspaper.result) # extract pe
kenya.newspaper.vc <- vcov(kenya.newspaper.result) # extract vcov
kenya.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
kenya.newspaper.scen <- cfName(kenya.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
kenya.newspaper.scen <- cfChange(kenya.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
kenya.newspaper.simbetas <- mvrnorm(sims, kenya.newspaper.pe, kenya.newspaper.vc) #create simbetas
kenya.newspaper.qoi <- logitsimfd(kenya.newspaper.scen, kenya.newspaper.simbetas, ci=0.95) #and get quantities of interest
kenya.newspaper.qoi.df <- as.data.frame(kenya.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "kenya")

#lesotho
mdata <- extractdata(newspaper.country, lesotho_data, na.rm = TRUE) #only extract data used in model
lesotho.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
lesotho.newspaper.pe <- fixef(lesotho.newspaper.result) # extract pe
lesotho.newspaper.vc <- vcov(lesotho.newspaper.result) # extract vcov
lesotho.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
lesotho.newspaper.scen <- cfName(lesotho.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
lesotho.newspaper.scen <- cfChange(lesotho.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
lesotho.newspaper.simbetas <- mvrnorm(sims, lesotho.newspaper.pe, lesotho.newspaper.vc) #create simbetas
lesotho.newspaper.qoi <- logitsimfd(lesotho.newspaper.scen, lesotho.newspaper.simbetas, ci=0.95) #and get quantities of interest
lesotho.newspaper.qoi.df <- as.data.frame(lesotho.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "lesotho")

#liberia
mdata <- extractdata(newspaper.country, liberia_data, na.rm = TRUE) #only extract data used in model
liberia.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
liberia.newspaper.pe <- fixef(liberia.newspaper.result) # extract pe
liberia.newspaper.vc <- vcov(liberia.newspaper.result) # extract vcov
liberia.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
liberia.newspaper.scen <- cfName(liberia.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
liberia.newspaper.scen <- cfChange(liberia.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
liberia.newspaper.simbetas <- mvrnorm(sims, liberia.newspaper.pe, liberia.newspaper.vc) #create simbetas
liberia.newspaper.qoi <- logitsimfd(liberia.newspaper.scen, liberia.newspaper.simbetas, ci=0.95) #and get quantities of interest
liberia.newspaper.qoi.df <- as.data.frame(liberia.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "liberia")

#madagascar
mdata <- extractdata(newspaper.country, madag_data, na.rm = TRUE) #only extract data used in model
madag.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
madag.newspaper.pe <- fixef(madag.newspaper.result) # extract pe
madag.newspaper.vc <- vcov(madag.newspaper.result) # extract vcov
madag.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
madag.newspaper.scen <- cfName(madag.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
madag.newspaper.scen <- cfChange(madag.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
madag.newspaper.simbetas <- mvrnorm(sims, madag.newspaper.pe, madag.newspaper.vc) #create simbetas
madag.newspaper.qoi <- logitsimfd(madag.newspaper.scen, madag.newspaper.simbetas, ci=0.95) #and get quantities of interest
madag.newspaper.qoi.df <- as.data.frame(madag.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "madagascar")

#malawi
mdata <- extractdata(newspaper.country, malawi_data, na.rm = TRUE) #only extract data used in model
malawi.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
malawi.newspaper.pe <- fixef(malawi.newspaper.result) # extract pe
malawi.newspaper.vc <- vcov(malawi.newspaper.result) # extract vcov
malawi.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
malawi.newspaper.scen <- cfName(malawi.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
malawi.newspaper.scen <- cfChange(malawi.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
malawi.newspaper.simbetas <- mvrnorm(sims, malawi.newspaper.pe, malawi.newspaper.vc) #create simbetas
malawi.newspaper.qoi <- logitsimfd(malawi.newspaper.scen, malawi.newspaper.simbetas, ci=0.95) #and get quantities of interest
malawi.newspaper.qoi.df <- as.data.frame(malawi.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "malawi")

#mali
mdata <- extractdata(newspaper.country, mali_data, na.rm = TRUE) #only extract data used in model
mali.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                        nAGQ = 0 ) # optimizer for faster but less precise fit
mali.newspaper.pe <- fixef(mali.newspaper.result) # extract pe
mali.newspaper.vc <- vcov(mali.newspaper.result) # extract vcov
mali.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mali.newspaper.scen <- cfName(mali.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
mali.newspaper.scen <- cfChange(mali.newspaper.scen, "newspaper", 
                         x = 5,
                         xpre=1,
                         scen=1)
mali.newspaper.simbetas <- mvrnorm(sims, mali.newspaper.pe, mali.newspaper.vc) #create simbetas
mali.newspaper.qoi <- logitsimfd(mali.newspaper.scen, mali.newspaper.simbetas, ci=0.95) #and get quantities of interest
mali.newspaper.qoi.df <- as.data.frame(mali.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "mali")

#mauritius
mdata <- extractdata(newspaper.country, maurit_data, na.rm = TRUE) #only extract data used in model
maurit.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
maurit.newspaper.pe <- fixef(maurit.newspaper.result) # extract pe
maurit.newspaper.vc <- vcov(maurit.newspaper.result) # extract vcov
maurit.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
maurit.newspaper.scen <- cfName(maurit.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
maurit.newspaper.scen <- cfChange(maurit.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
maurit.newspaper.simbetas <- mvrnorm(sims, maurit.newspaper.pe, maurit.newspaper.vc) #create simbetas
maurit.newspaper.qoi <- logitsimfd(maurit.newspaper.scen, maurit.newspaper.simbetas, ci=0.95) #and get quantities of interest
maurit.newspaper.qoi.df <- as.data.frame(maurit.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "mauritius")

#morocco
mdata <- extractdata(newspaper.country, morocco_data, na.rm = TRUE) #only extract data used in model
morocco.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
morocco.newspaper.pe <- fixef(morocco.newspaper.result) # extract pe
morocco.newspaper.vc <- vcov(morocco.newspaper.result) # extract vcov
morocco.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
morocco.newspaper.scen <- cfName(morocco.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
morocco.newspaper.scen <- cfChange(morocco.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
morocco.newspaper.simbetas <- mvrnorm(sims, morocco.newspaper.pe, morocco.newspaper.vc) #create simbetas
morocco.newspaper.qoi <- logitsimfd(morocco.newspaper.scen, morocco.newspaper.simbetas, ci=0.95) #and get quantities of interest
morocco.newspaper.qoi.df <- as.data.frame(morocco.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "morocco")

#mozambique
mdata <- extractdata(newspaper.country, mozamb_data, na.rm = TRUE) #only extract data used in model
mozamb.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
mozamb.newspaper.pe <- fixef(mozamb.newspaper.result) # extract pe
mozamb.newspaper.vc <- vcov(mozamb.newspaper.result) # extract vcov
mozamb.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mozamb.newspaper.scen <- cfName(mozamb.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
mozamb.newspaper.scen <- cfChange(mozamb.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
mozamb.newspaper.simbetas <- mvrnorm(sims, mozamb.newspaper.pe, mozamb.newspaper.vc) #create simbetas
mozamb.newspaper.qoi <- logitsimfd(mozamb.newspaper.scen, mozamb.newspaper.simbetas, ci=0.95) #and get quantities of interest
mozamb.newspaper.qoi.df <- as.data.frame(mozamb.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "mozambique")

#namibia
mdata <- extractdata(newspaper.country, namibia_data, na.rm = TRUE) #only extract data used in model
namibia.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
namibia.newspaper.pe <- fixef(namibia.newspaper.result) # extract pe
namibia.newspaper.vc <- vcov(namibia.newspaper.result) # extract vcov
namibia.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
namibia.newspaper.scen <- cfName(namibia.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
namibia.newspaper.scen <- cfChange(namibia.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
namibia.newspaper.simbetas <- mvrnorm(sims, namibia.newspaper.pe, namibia.newspaper.vc) #create simbetas
namibia.newspaper.qoi <- logitsimfd(namibia.newspaper.scen, namibia.newspaper.simbetas, ci=0.95) #and get quantities of interest
namibia.newspaper.qoi.df <- as.data.frame(namibia.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "namibia")

#niger
mdata <- extractdata(newspaper.country, niger_data, na.rm = TRUE) #only extract data used in model
niger.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                         nAGQ = 0 ) # optimizer for faster but less precise fit
niger.newspaper.pe <- fixef(niger.newspaper.result) # extract pe
niger.newspaper.vc <- vcov(niger.newspaper.result) # extract vcov
niger.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
niger.newspaper.scen <- cfName(niger.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
niger.newspaper.scen <- cfChange(niger.newspaper.scen, "newspaper", 
                          x = 5,
                          xpre=1,
                          scen=1)
niger.newspaper.simbetas <- mvrnorm(sims, niger.newspaper.pe, niger.newspaper.vc) #create simbetas
niger.newspaper.qoi <- logitsimfd(niger.newspaper.scen, niger.newspaper.simbetas, ci=0.95) #and get quantities of interest
niger.newspaper.qoi.df <- as.data.frame(niger.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "niger")

#nigeria
mdata <- extractdata(newspaper.country, nigeria_data, na.rm = TRUE) #only extract data used in model
nigeria.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
nigeria.newspaper.pe <- fixef(nigeria.newspaper.result) # extract pe
nigeria.newspaper.vc <- vcov(nigeria.newspaper.result) # extract vcov
nigeria.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
nigeria.newspaper.scen <- cfName(nigeria.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
nigeria.newspaper.scen <- cfChange(nigeria.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
nigeria.newspaper.simbetas <- mvrnorm(sims, nigeria.newspaper.pe, nigeria.newspaper.vc) #create simbetas
nigeria.newspaper.qoi <- logitsimfd(nigeria.newspaper.scen, nigeria.newspaper.simbetas, ci=0.95) #and get quantities of interest
nigeria.newspaper.qoi.df <- as.data.frame(nigeria.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "nigeria")

#são tomé and príncipe
mdata <- extractdata(newspaper.country, stp_data, na.rm = TRUE) #only extract data used in model
stp.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                       nAGQ = 0 ) # optimizer for faster but less precise fit
stp.newspaper.pe <- fixef(stp.newspaper.result) # extract pe
stp.newspaper.vc <- vcov(stp.newspaper.result) # extract vcov
stp.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
stp.newspaper.scen <- cfName(stp.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
stp.newspaper.scen <- cfChange(stp.newspaper.scen, "newspaper", 
                        x = 5,
                        xpre=1,
                        scen=1)
stp.newspaper.simbetas <- mvrnorm(sims, stp.newspaper.pe, stp.newspaper.vc) #create simbetas
stp.newspaper.qoi <- logitsimfd(stp.newspaper.scen, stp.newspaper.simbetas, ci=0.95) #and get quantities of interest
stp.newspaper.qoi.df <- as.data.frame(stp.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "são tomé and príncipe")

#senegal
mdata <- extractdata(newspaper.country, senegal_data, na.rm = TRUE) #only extract data used in model
senegal.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
senegal.newspaper.pe <- fixef(senegal.newspaper.result) # extract pe
senegal.newspaper.vc <- vcov(senegal.newspaper.result) # extract vcov
senegal.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
senegal.newspaper.scen <- cfName(senegal.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
senegal.newspaper.scen <- cfChange(senegal.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
senegal.newspaper.simbetas <- mvrnorm(sims, senegal.newspaper.pe, senegal.newspaper.vc) #create simbetas
senegal.newspaper.qoi <- logitsimfd(senegal.newspaper.scen, senegal.newspaper.simbetas, ci=0.95) #and get quantities of interest
senegal.newspaper.qoi.df <- as.data.frame(senegal.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "senegal")

#sierra
mdata <- extractdata(newspaper.country, sierra_data, na.rm = TRUE) #only extract data used in model
sierra.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
sierra.newspaper.pe <- fixef(sierra.newspaper.result) # extract pe
sierra.newspaper.vc <- vcov(sierra.newspaper.result) # extract vcov
sierra.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
sierra.newspaper.scen <- cfName(sierra.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
sierra.newspaper.scen <- cfChange(sierra.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
sierra.newspaper.simbetas <- mvrnorm(sims, sierra.newspaper.pe, sierra.newspaper.vc) #create simbetas
sierra.newspaper.qoi <- logitsimfd(sierra.newspaper.scen, sierra.newspaper.simbetas, ci=0.95) #and get quantities of interest
sierra.newspaper.qoi.df <- as.data.frame(sierra.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "sierra leone")

#southaf
mdata <- extractdata(newspaper.country, southaf_data, na.rm = TRUE) #only extract data used in model
southaf.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
southaf.newspaper.pe <- fixef(southaf.newspaper.result) # extract pe
southaf.newspaper.vc <- vcov(southaf.newspaper.result) # extract vcov
southaf.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
southaf.newspaper.scen <- cfName(southaf.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
southaf.newspaper.scen <- cfChange(southaf.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
southaf.newspaper.simbetas <- mvrnorm(sims, southaf.newspaper.pe, southaf.newspaper.vc) #create simbetas
southaf.newspaper.qoi <- logitsimfd(southaf.newspaper.scen, southaf.newspaper.simbetas, ci=0.95) #and get quantities of interest
southaf.newspaper.qoi.df <- as.data.frame(southaf.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "south africa")

#swaz
mdata <- extractdata(newspaper.country, swaz_data, na.rm = TRUE) #only extract data used in model
swaz.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                        nAGQ = 0 ) # optimizer for faster but less precise fit
swaz.newspaper.pe <- fixef(swaz.newspaper.result) # extract pe
swaz.newspaper.vc <- vcov(swaz.newspaper.result) # extract vcov
swaz.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
swaz.newspaper.scen <- cfName(swaz.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
swaz.newspaper.scen <- cfChange(swaz.newspaper.scen, "newspaper", 
                         x = 5,
                         xpre=1,
                         scen=1)
swaz.newspaper.simbetas <- mvrnorm(sims, swaz.newspaper.pe, swaz.newspaper.vc) #create simbetas
swaz.newspaper.qoi <- logitsimfd(swaz.newspaper.scen, swaz.newspaper.simbetas, ci=0.95) #and get quantities of interest
swaz.newspaper.qoi.df <- as.data.frame(swaz.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "swaziland")

#tanz
mdata <- extractdata(newspaper.country, tanz_data, na.rm = TRUE) #only extract data used in model
tanz.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                        nAGQ = 0 ) # optimizer for faster but less precise fit
tanz.newspaper.pe <- fixef(tanz.newspaper.result) # extract pe
tanz.newspaper.vc <- vcov(tanz.newspaper.result) # extract vcov
tanz.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tanz.newspaper.scen <- cfName(tanz.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
tanz.newspaper.scen <- cfChange(tanz.newspaper.scen, "newspaper", 
                         x = 5,
                         xpre=1,
                         scen=1)
tanz.newspaper.simbetas <- mvrnorm(sims, tanz.newspaper.pe, tanz.newspaper.vc) #create simbetas
tanz.newspaper.qoi <- logitsimfd(tanz.newspaper.scen, tanz.newspaper.simbetas, ci=0.95) #and get quantities of interest
tanz.newspaper.qoi.df <- as.data.frame(tanz.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "tanzania")

#togo
mdata <- extractdata(newspaper.country, togo_data, na.rm = TRUE) #only extract data used in model
togo.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                        nAGQ = 0 ) # optimizer for faster but less precise fit
togo.newspaper.pe <- fixef(togo.newspaper.result) # extract pe
togo.newspaper.vc <- vcov(togo.newspaper.result) # extract vcov
togo.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
togo.newspaper.scen <- cfName(togo.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
togo.newspaper.scen <- cfChange(togo.newspaper.scen, "newspaper", 
                         x = 5,
                         xpre=1,
                         scen=1)
togo.newspaper.simbetas <- mvrnorm(sims, togo.newspaper.pe, togo.newspaper.vc) #create simbetas
togo.newspaper.qoi <- logitsimfd(togo.newspaper.scen, togo.newspaper.simbetas, ci=0.95) #and get quantities of interest
togo.newspaper.qoi.df <- as.data.frame(togo.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "togo")

#tunisia
mdata <- extractdata(newspaper.country, tunisia_data, na.rm = TRUE) #only extract data used in model
tunisia.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                           nAGQ = 0 ) # optimizer for faster but less precise fit
tunisia.newspaper.pe <- fixef(tunisia.newspaper.result) # extract pe
tunisia.newspaper.vc <- vcov(tunisia.newspaper.result) # extract vcov
tunisia.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tunisia.newspaper.scen <- cfName(tunisia.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
tunisia.newspaper.scen <- cfChange(tunisia.newspaper.scen, "newspaper", 
                            x = 5,
                            xpre=1,
                            scen=1)
tunisia.newspaper.simbetas <- mvrnorm(sims, tunisia.newspaper.pe, tunisia.newspaper.vc) #create simbetas
tunisia.newspaper.qoi <- logitsimfd(tunisia.newspaper.scen, tunisia.newspaper.simbetas, ci=0.95) #and get quantities of interest
tunisia.newspaper.qoi.df <- as.data.frame(tunisia.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "tunisia")

#uganda
mdata <- extractdata(newspaper.country, uganda_data, na.rm = TRUE) #only extract data used in model
uganda.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
uganda.newspaper.pe <- fixef(uganda.newspaper.result) # extract pe
uganda.newspaper.vc <- vcov(uganda.newspaper.result) # extract vcov
uganda.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
uganda.newspaper.scen <- cfName(uganda.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
uganda.newspaper.scen <- cfChange(uganda.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
uganda.newspaper.simbetas <- mvrnorm(sims, uganda.newspaper.pe, uganda.newspaper.vc) #create simbetas
uganda.newspaper.qoi <- logitsimfd(uganda.newspaper.scen, uganda.newspaper.simbetas, ci=0.95) #and get quantities of interest
uganda.newspaper.qoi.df <- as.data.frame(uganda.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "uganda")

#zambia
mdata <- extractdata(newspaper.country, zambia_data, na.rm = TRUE) #only extract data used in model
zambia.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
zambia.newspaper.pe <- fixef(zambia.newspaper.result) # extract pe
zambia.newspaper.vc <- vcov(zambia.newspaper.result) # extract vcov
zambia.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zambia.newspaper.scen <- cfName(zambia.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
zambia.newspaper.scen <- cfChange(zambia.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
zambia.newspaper.simbetas <- mvrnorm(sims, zambia.newspaper.pe, zambia.newspaper.vc) #create simbetas
zambia.newspaper.qoi <- logitsimfd(zambia.newspaper.scen, zambia.newspaper.simbetas, ci=0.95) #and get quantities of interest
zambia.newspaper.qoi.df <- as.data.frame(zambia.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "zambia")

#zimbab
mdata <- extractdata(newspaper.country, zimbab_data, na.rm = TRUE) #only extract data used in model
zimbab.newspaper.result <- glmer(newspaper.country, family=binomial, data=mdata, #estimate the model
                          nAGQ = 0 ) # optimizer for faster but less precise fit
zimbab.newspaper.pe <- fixef(zimbab.newspaper.result) # extract pe
zimbab.newspaper.vc <- vcov(zimbab.newspaper.result) # extract vcov
zimbab.newspaper.scen <- cfMake(newspaper.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zimbab.newspaper.scen <- cfName(zimbab.newspaper.scen, "newspaper 1 vs 5", scen=1) #configure scenario of 1 to 5
zimbab.newspaper.scen <- cfChange(zimbab.newspaper.scen, "newspaper", 
                           x = 5,
                           xpre=1,
                           scen=1)
zimbab.newspaper.simbetas <- mvrnorm(sims, zimbab.newspaper.pe, zimbab.newspaper.vc) #create simbetas
zimbab.newspaper.qoi <- logitsimfd(zimbab.newspaper.scen, zimbab.newspaper.simbetas, ci=0.95) #and get quantities of interest
zimbab.newspaper.qoi.df <- as.data.frame(zimbab.newspaper.qoi) %>% mutate(medium = 'newspaper', country = "zimbabwe")


###########
# Internet
###########

internet.country <- #use this model to estimate result
  sexuality2 ~
  internet +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 | district)  #add VI for District

internet.country.nodist <- #use this mdoel to create scenarios (drop district)
  sexuality2 ~
  internet +
  media_nointernet + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 

#Benin
mdata <- extractdata(internet.country, benin_data, na.rm = TRUE) #only extract data used in model
benin.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
benin.internet.pe <- fixef(benin.internet.result) # extract pe
benin.internet.vc <- vcov(benin.internet.result) # extract vcov
benin.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
benin.internet.scen <- cfName(benin.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
benin.internet.scen <- cfChange(benin.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
benin.internet.simbetas <- mvrnorm(sims, benin.internet.pe, benin.internet.vc) #create simbetas
benin.internet.qoi <- logitsimfd(benin.internet.scen, benin.internet.simbetas, ci=0.95) #and get quantities of interest
benin.internet.qoi.df <- as.data.frame(benin.internet.qoi) %>% mutate(medium = 'internet', country = "benin") #create df from qoi

#Botswana
mdata <- extractdata(internet.country.nodist, botswana_data, na.rm = TRUE) #use nodist model for botswana b/c there are no districts (i.e. no random effects)
botswana.internet.result <- glm(internet.country.nodist, family=binomial, data=mdata) #estimate the model
botswana.internet.pe <-botswana.internet.result$coefficients  # extract pe
botswana.internet.vc <- vcov(botswana.internet.result) # extract vcov
botswana.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
botswana.internet.scen <- cfName(botswana.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
botswana.internet.scen <- cfChange(botswana.internet.scen, "internet", 
                                x = 5,
                                xpre=1,
                                scen=1)
botswana.internet.simbetas <- mvrnorm(sims, botswana.internet.pe, botswana.internet.vc) #create simbetas
botswana.internet.qoi <- logitsimfd(botswana.internet.scen, botswana.internet.simbetas, ci=0.95) #and get quantities of interest
botswana.internet.qoi.df <- as.data.frame(botswana.internet.qoi) %>% mutate(medium = 'internet', country = "botswana")

#burkina faso
mdata <- extractdata(internet.country, burkina_data, na.rm = TRUE) #only extract data used in model
burkina.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
burkina.internet.pe <- fixef(burkina.internet.result) # extract pe
burkina.internet.vc <- vcov(burkina.internet.result) # extract vcov
burkina.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burkina.internet.scen <- cfName(burkina.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
burkina.internet.scen <- cfChange(burkina.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
burkina.internet.simbetas <- mvrnorm(sims, burkina.internet.pe, burkina.internet.vc) #create simbetas
burkina.internet.qoi <- logitsimfd(burkina.internet.scen, burkina.internet.simbetas, ci=0.95) #and get quantities of interest
burkina.internet.qoi.df <- as.data.frame(burkina.internet.qoi) %>% mutate(medium = 'internet', country = "burkina faso")

#burundi
mdata <- extractdata(internet.country, burundi_data, na.rm = TRUE) #only extract data used in model
burundi.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
burundi.internet.pe <- fixef(burundi.internet.result) # extract pe
burundi.internet.vc <- vcov(burundi.internet.result) # extract vcov
burundi.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burundi.internet.scen <- cfName(burundi.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
burundi.internet.scen <- cfChange(burundi.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
burundi.internet.simbetas <- mvrnorm(sims, burundi.internet.pe, burundi.internet.vc) #create simbetas
burundi.internet.qoi <- logitsimfd(burundi.internet.scen, burundi.internet.simbetas, ci=0.95) #and get quantities of interest
burundi.internet.qoi.df <- as.data.frame(burundi.internet.qoi) %>% mutate(medium = 'internet', country = "burundi")

#cameroon
mdata <- extractdata(internet.country, cameroon_data, na.rm = TRUE) #only extract data used in model
cameroon.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                   nAGQ = 0 ) # optimizer for faster but less precise fit
cameroon.internet.pe <- fixef(cameroon.internet.result) # extract pe
cameroon.internet.vc <- vcov(cameroon.internet.result) # extract vcov
cameroon.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cameroon.internet.scen <- cfName(cameroon.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
cameroon.internet.scen <- cfChange(cameroon.internet.scen, "internet", 
                                    x = 5,
                                    xpre=1,
                                    scen=1)
cameroon.internet.simbetas <- mvrnorm(sims, cameroon.internet.pe, cameroon.internet.vc) #create simbetas
cameroon.internet.qoi <- logitsimfd(cameroon.internet.scen, cameroon.internet.simbetas, ci=0.95) #and get quantities of interest
cameroon.internet.qoi.df <- as.data.frame(cameroon.internet.qoi) %>% mutate(medium = 'internet', country = "cameroon")

#cape verde
mdata <- extractdata(internet.country, capeverde_data, na.rm = TRUE) #only extract data used in model
capeverde.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                    nAGQ = 0 ) # optimizer for faster but less precise fit
capeverde.internet.pe <- fixef(capeverde.internet.result) # extract pe
capeverde.internet.vc <- vcov(capeverde.internet.result) # extract vcov
capeverde.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
capeverde.internet.scen <- cfName(capeverde.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
capeverde.internet.scen <- cfChange(capeverde.internet.scen, "internet", 
                                     x = 5,
                                     xpre=1,
                                     scen=1)
capeverde.internet.simbetas <- mvrnorm(sims, capeverde.internet.pe, capeverde.internet.vc) #create simbetas
capeverde.internet.qoi <- logitsimfd(capeverde.internet.scen, capeverde.internet.simbetas, ci=0.95) #and get quantities of interest
capeverde.internet.qoi.df <- as.data.frame(capeverde.internet.qoi) %>% mutate(medium = 'internet', country = "cape verde")

#cote d'ivoire
mdata <- extractdata(internet.country, cdi_data, na.rm = TRUE) #only extract data used in model
cdi.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
cdi.internet.pe <- fixef(cdi.internet.result) # extract pe
cdi.internet.vc <- vcov(cdi.internet.result) # extract vcov
cdi.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cdi.internet.scen <- cfName(cdi.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
cdi.internet.scen <- cfChange(cdi.internet.scen, "internet", 
                               x = 5,
                               xpre=1,
                               scen=1)
cdi.internet.simbetas <- mvrnorm(sims, cdi.internet.pe, cdi.internet.vc) #create simbetas
cdi.internet.qoi <- logitsimfd(cdi.internet.scen, cdi.internet.simbetas, ci=0.95) #and get quantities of interest
cdi.internet.qoi.df <- as.data.frame(cdi.internet.qoi) %>% mutate(medium = 'internet', country = "cote d'ivoire")

#gabon
mdata <- extractdata(internet.country, gabon_data, na.rm = TRUE) #only extract data used in model
gabon.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
gabon.internet.pe <- fixef(gabon.internet.result) # extract pe
gabon.internet.vc <- vcov(gabon.internet.result) # extract vcov
gabon.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
gabon.internet.scen <- cfName(gabon.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
gabon.internet.scen <- cfChange(gabon.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
gabon.internet.simbetas <- mvrnorm(sims, gabon.internet.pe, gabon.internet.vc) #create simbetas
gabon.internet.qoi <- logitsimfd(gabon.internet.scen, gabon.internet.simbetas, ci=0.95) #and get quantities of interest
gabon.internet.qoi.df <- as.data.frame(gabon.internet.qoi) %>% mutate(medium = 'internet', country = "gabon")

#ghana
mdata <- extractdata(internet.country, ghana_data, na.rm = TRUE) #only extract data used in model
ghana.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
ghana.internet.pe <- fixef(ghana.internet.result) # extract pe
ghana.internet.vc <- vcov(ghana.internet.result) # extract vcov
ghana.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
ghana.internet.scen <- cfName(ghana.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
ghana.internet.scen <- cfChange(ghana.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
ghana.internet.simbetas <- mvrnorm(sims, ghana.internet.pe, ghana.internet.vc) #create simbetas
ghana.internet.qoi <- logitsimfd(ghana.internet.scen, ghana.internet.simbetas, ci=0.95) #and get quantities of interest
ghana.internet.qoi.df <- as.data.frame(ghana.internet.qoi) %>% mutate(medium = 'internet', country = "ghana")

#guinea 
mdata <- extractdata(internet.country, guinea_data, na.rm = TRUE) #only extract data used in model
guinea.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
guinea.internet.pe <- fixef(guinea.internet.result) # extract pe
guinea.internet.vc <- vcov(guinea.internet.result) # extract vcov
guinea.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
guinea.internet.scen <- cfName(guinea.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
guinea.internet.scen <- cfChange(guinea.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
guinea.internet.simbetas <- mvrnorm(sims, guinea.internet.pe, guinea.internet.vc) #create simbetas
guinea.internet.qoi <- logitsimfd(guinea.internet.scen, guinea.internet.simbetas, ci=0.95) #and get quantities of interest
guinea.internet.qoi.df <- as.data.frame(guinea.internet.qoi) %>% mutate(medium = 'internet', country = "guinea")

#kenya
mdata <- extractdata(internet.country, kenya_data, na.rm = TRUE) #only extract data used in model
kenya.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
kenya.internet.pe <- fixef(kenya.internet.result) # extract pe
kenya.internet.vc <- vcov(kenya.internet.result) # extract vcov
kenya.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
kenya.internet.scen <- cfName(kenya.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
kenya.internet.scen <- cfChange(kenya.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
kenya.internet.simbetas <- mvrnorm(sims, kenya.internet.pe, kenya.internet.vc) #create simbetas
kenya.internet.qoi <- logitsimfd(kenya.internet.scen, kenya.internet.simbetas, ci=0.95) #and get quantities of interest
kenya.internet.qoi.df <- as.data.frame(kenya.internet.qoi) %>% mutate(medium = 'internet', country = "kenya")

#lesotho
mdata <- extractdata(internet.country, lesotho_data, na.rm = TRUE) #only extract data used in model
lesotho.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
lesotho.internet.pe <- fixef(lesotho.internet.result) # extract pe
lesotho.internet.vc <- vcov(lesotho.internet.result) # extract vcov
lesotho.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
lesotho.internet.scen <- cfName(lesotho.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
lesotho.internet.scen <- cfChange(lesotho.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
lesotho.internet.simbetas <- mvrnorm(sims, lesotho.internet.pe, lesotho.internet.vc) #create simbetas
lesotho.internet.qoi <- logitsimfd(lesotho.internet.scen, lesotho.internet.simbetas, ci=0.95) #and get quantities of interest
lesotho.internet.qoi.df <- as.data.frame(lesotho.internet.qoi) %>% mutate(medium = 'internet', country = "lesotho")

#liberia
mdata <- extractdata(internet.country, liberia_data, na.rm = TRUE) #only extract data used in model
liberia.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
liberia.internet.pe <- fixef(liberia.internet.result) # extract pe
liberia.internet.vc <- vcov(liberia.internet.result) # extract vcov
liberia.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
liberia.internet.scen <- cfName(liberia.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
liberia.internet.scen <- cfChange(liberia.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
liberia.internet.simbetas <- mvrnorm(sims, liberia.internet.pe, liberia.internet.vc) #create simbetas
liberia.internet.qoi <- logitsimfd(liberia.internet.scen, liberia.internet.simbetas, ci=0.95) #and get quantities of interest
liberia.internet.qoi.df <- as.data.frame(liberia.internet.qoi) %>% mutate(medium = 'internet', country = "liberia")

#madagascar
mdata <- extractdata(internet.country, madag_data, na.rm = TRUE) #only extract data used in model
madag.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
madag.internet.pe <- fixef(madag.internet.result) # extract pe
madag.internet.vc <- vcov(madag.internet.result) # extract vcov
madag.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
madag.internet.scen <- cfName(madag.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
madag.internet.scen <- cfChange(madag.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
madag.internet.simbetas <- mvrnorm(sims, madag.internet.pe, madag.internet.vc) #create simbetas
madag.internet.qoi <- logitsimfd(madag.internet.scen, madag.internet.simbetas, ci=0.95) #and get quantities of interest
madag.internet.qoi.df <- as.data.frame(madag.internet.qoi) %>% mutate(medium = 'internet', country = "madagascar")

#malawi
mdata <- extractdata(internet.country, malawi_data, na.rm = TRUE) #only extract data used in model
malawi.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
malawi.internet.pe <- fixef(malawi.internet.result) # extract pe
malawi.internet.vc <- vcov(malawi.internet.result) # extract vcov
malawi.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
malawi.internet.scen <- cfName(malawi.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
malawi.internet.scen <- cfChange(malawi.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
malawi.internet.simbetas <- mvrnorm(sims, malawi.internet.pe, malawi.internet.vc) #create simbetas
malawi.internet.qoi <- logitsimfd(malawi.internet.scen, malawi.internet.simbetas, ci=0.95) #and get quantities of interest
malawi.internet.qoi.df <- as.data.frame(malawi.internet.qoi) %>% mutate(medium = 'internet', country = "malawi")

#mali
mdata <- extractdata(internet.country, mali_data, na.rm = TRUE) #only extract data used in model
mali.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
mali.internet.pe <- fixef(mali.internet.result) # extract pe
mali.internet.vc <- vcov(mali.internet.result) # extract vcov
mali.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mali.internet.scen <- cfName(mali.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
mali.internet.scen <- cfChange(mali.internet.scen, "internet", 
                                x = 5,
                                xpre=1,
                                scen=1)
mali.internet.simbetas <- mvrnorm(sims, mali.internet.pe, mali.internet.vc) #create simbetas
mali.internet.qoi <- logitsimfd(mali.internet.scen, mali.internet.simbetas, ci=0.95) #and get quantities of interest
mali.internet.qoi.df <- as.data.frame(mali.internet.qoi) %>% mutate(medium = 'internet', country = "mali")

#mauritius
mdata <- extractdata(internet.country, maurit_data, na.rm = TRUE) #only extract data used in model
maurit.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
maurit.internet.pe <- fixef(maurit.internet.result) # extract pe
maurit.internet.vc <- vcov(maurit.internet.result) # extract vcov
maurit.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
maurit.internet.scen <- cfName(maurit.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
maurit.internet.scen <- cfChange(maurit.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
maurit.internet.simbetas <- mvrnorm(sims, maurit.internet.pe, maurit.internet.vc) #create simbetas
maurit.internet.qoi <- logitsimfd(maurit.internet.scen, maurit.internet.simbetas, ci=0.95) #and get quantities of interest
maurit.internet.qoi.df <- as.data.frame(maurit.internet.qoi) %>% mutate(medium = 'internet', country = "mauritius")

#morocco
mdata <- extractdata(internet.country, morocco_data, na.rm = TRUE) #only extract data used in model
morocco.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
morocco.internet.pe <- fixef(morocco.internet.result) # extract pe
morocco.internet.vc <- vcov(morocco.internet.result) # extract vcov
morocco.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
morocco.internet.scen <- cfName(morocco.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
morocco.internet.scen <- cfChange(morocco.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
morocco.internet.simbetas <- mvrnorm(sims, morocco.internet.pe, morocco.internet.vc) #create simbetas
morocco.internet.qoi <- logitsimfd(morocco.internet.scen, morocco.internet.simbetas, ci=0.95) #and get quantities of interest
morocco.internet.qoi.df <- as.data.frame(morocco.internet.qoi) %>% mutate(medium = 'internet', country = "morocco")

#mozambique
mdata <- extractdata(internet.country, mozamb_data, na.rm = TRUE) #only extract data used in model
mozamb.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
mozamb.internet.pe <- fixef(mozamb.internet.result) # extract pe
mozamb.internet.vc <- vcov(mozamb.internet.result) # extract vcov
mozamb.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mozamb.internet.scen <- cfName(mozamb.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
mozamb.internet.scen <- cfChange(mozamb.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
mozamb.internet.simbetas <- mvrnorm(sims, mozamb.internet.pe, mozamb.internet.vc) #create simbetas
mozamb.internet.qoi <- logitsimfd(mozamb.internet.scen, mozamb.internet.simbetas, ci=0.95) #and get quantities of interest
mozamb.internet.qoi.df <- as.data.frame(mozamb.internet.qoi) %>% mutate(medium = 'internet', country = "mozambique")

#namibia
mdata <- extractdata(internet.country, namibia_data, na.rm = TRUE) #only extract data used in model
namibia.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
namibia.internet.pe <- fixef(namibia.internet.result) # extract pe
namibia.internet.vc <- vcov(namibia.internet.result) # extract vcov
namibia.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
namibia.internet.scen <- cfName(namibia.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
namibia.internet.scen <- cfChange(namibia.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
namibia.internet.simbetas <- mvrnorm(sims, namibia.internet.pe, namibia.internet.vc) #create simbetas
namibia.internet.qoi <- logitsimfd(namibia.internet.scen, namibia.internet.simbetas, ci=0.95) #and get quantities of interest
namibia.internet.qoi.df <- as.data.frame(namibia.internet.qoi) %>% mutate(medium = 'internet', country = "namibia")

#niger
mdata <- extractdata(internet.country, niger_data, na.rm = TRUE) #only extract data used in model
niger.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
niger.internet.pe <- fixef(niger.internet.result) # extract pe
niger.internet.vc <- vcov(niger.internet.result) # extract vcov
niger.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
niger.internet.scen <- cfName(niger.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
niger.internet.scen <- cfChange(niger.internet.scen, "internet", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
niger.internet.simbetas <- mvrnorm(sims, niger.internet.pe, niger.internet.vc) #create simbetas
niger.internet.qoi <- logitsimfd(niger.internet.scen, niger.internet.simbetas, ci=0.95) #and get quantities of interest
niger.internet.qoi.df <- as.data.frame(niger.internet.qoi) %>% mutate(medium = 'internet', country = "niger")

#nigeria
mdata <- extractdata(internet.country, nigeria_data, na.rm = TRUE) #only extract data used in model
nigeria.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
nigeria.internet.pe <- fixef(nigeria.internet.result) # extract pe
nigeria.internet.vc <- vcov(nigeria.internet.result) # extract vcov
nigeria.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
nigeria.internet.scen <- cfName(nigeria.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
nigeria.internet.scen <- cfChange(nigeria.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
nigeria.internet.simbetas <- mvrnorm(sims, nigeria.internet.pe, nigeria.internet.vc) #create simbetas
nigeria.internet.qoi <- logitsimfd(nigeria.internet.scen, nigeria.internet.simbetas, ci=0.95) #and get quantities of interest
nigeria.internet.qoi.df <- as.data.frame(nigeria.internet.qoi) %>% mutate(medium = 'internet', country = "nigeria")

#são tomé and príncipe
mdata <- extractdata(internet.country, stp_data, na.rm = TRUE) #only extract data used in model
stp.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
stp.internet.pe <- fixef(stp.internet.result) # extract pe
stp.internet.vc <- vcov(stp.internet.result) # extract vcov
stp.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
stp.internet.scen <- cfName(stp.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
stp.internet.scen <- cfChange(stp.internet.scen, "internet", 
                               x = 5,
                               xpre=1,
                               scen=1)
stp.internet.simbetas <- mvrnorm(sims, stp.internet.pe, stp.internet.vc) #create simbetas
stp.internet.qoi <- logitsimfd(stp.internet.scen, stp.internet.simbetas, ci=0.95) #and get quantities of interest
stp.internet.qoi.df <- as.data.frame(stp.internet.qoi) %>% mutate(medium = 'internet', country = "são tomé and príncipe")

#senegal
mdata <- extractdata(internet.country, senegal_data, na.rm = TRUE) #only extract data used in model
senegal.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
senegal.internet.pe <- fixef(senegal.internet.result) # extract pe
senegal.internet.vc <- vcov(senegal.internet.result) # extract vcov
senegal.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
senegal.internet.scen <- cfName(senegal.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
senegal.internet.scen <- cfChange(senegal.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
senegal.internet.simbetas <- mvrnorm(sims, senegal.internet.pe, senegal.internet.vc) #create simbetas
senegal.internet.qoi <- logitsimfd(senegal.internet.scen, senegal.internet.simbetas, ci=0.95) #and get quantities of interest
senegal.internet.qoi.df <- as.data.frame(senegal.internet.qoi) %>% mutate(medium = 'internet', country = "senegal")

#sierra
mdata <- extractdata(internet.country, sierra_data, na.rm = TRUE) #only extract data used in model
sierra.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
sierra.internet.pe <- fixef(sierra.internet.result) # extract pe
sierra.internet.vc <- vcov(sierra.internet.result) # extract vcov
sierra.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
sierra.internet.scen <- cfName(sierra.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
sierra.internet.scen <- cfChange(sierra.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
sierra.internet.simbetas <- mvrnorm(sims, sierra.internet.pe, sierra.internet.vc) #create simbetas
sierra.internet.qoi <- logitsimfd(sierra.internet.scen, sierra.internet.simbetas, ci=0.95) #and get quantities of interest
sierra.internet.qoi.df <- as.data.frame(sierra.internet.qoi) %>% mutate(medium = 'internet', country = "sierra leone")

#southaf
mdata <- extractdata(internet.country, southaf_data, na.rm = TRUE) #only extract data used in model
southaf.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
southaf.internet.pe <- fixef(southaf.internet.result) # extract pe
southaf.internet.vc <- vcov(southaf.internet.result) # extract vcov
southaf.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
southaf.internet.scen <- cfName(southaf.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
southaf.internet.scen <- cfChange(southaf.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
southaf.internet.simbetas <- mvrnorm(sims, southaf.internet.pe, southaf.internet.vc) #create simbetas
southaf.internet.qoi <- logitsimfd(southaf.internet.scen, southaf.internet.simbetas, ci=0.95) #and get quantities of interest
southaf.internet.qoi.df <- as.data.frame(southaf.internet.qoi) %>% mutate(medium = 'internet', country = "south africa")

#swaz
mdata <- extractdata(internet.country, swaz_data, na.rm = TRUE) #only extract data used in model
swaz.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
swaz.internet.pe <- fixef(swaz.internet.result) # extract pe
swaz.internet.vc <- vcov(swaz.internet.result) # extract vcov
swaz.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
swaz.internet.scen <- cfName(swaz.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
swaz.internet.scen <- cfChange(swaz.internet.scen, "internet", 
                                x = 5,
                                xpre=1,
                                scen=1)
swaz.internet.simbetas <- mvrnorm(sims, swaz.internet.pe, swaz.internet.vc) #create simbetas
swaz.internet.qoi <- logitsimfd(swaz.internet.scen, swaz.internet.simbetas, ci=0.95) #and get quantities of interest
swaz.internet.qoi.df <- as.data.frame(swaz.internet.qoi) %>% mutate(medium = 'internet', country = "swaziland")

#tanz
mdata <- extractdata(internet.country, tanz_data, na.rm = TRUE) #only extract data used in model
tanz.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
tanz.internet.pe <- fixef(tanz.internet.result) # extract pe
tanz.internet.vc <- vcov(tanz.internet.result) # extract vcov
tanz.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tanz.internet.scen <- cfName(tanz.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
tanz.internet.scen <- cfChange(tanz.internet.scen, "internet", 
                                x = 5,
                                xpre=1,
                                scen=1)
tanz.internet.simbetas <- mvrnorm(sims, tanz.internet.pe, tanz.internet.vc) #create simbetas
tanz.internet.qoi <- logitsimfd(tanz.internet.scen, tanz.internet.simbetas, ci=0.95) #and get quantities of interest
tanz.internet.qoi.df <- as.data.frame(tanz.internet.qoi) %>% mutate(medium = 'internet', country = "tanzania")

#togo
mdata <- extractdata(internet.country, togo_data, na.rm = TRUE) #only extract data used in model
togo.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
togo.internet.pe <- fixef(togo.internet.result) # extract pe
togo.internet.vc <- vcov(togo.internet.result) # extract vcov
togo.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
togo.internet.scen <- cfName(togo.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
togo.internet.scen <- cfChange(togo.internet.scen, "internet", 
                                x = 5,
                                xpre=1,
                                scen=1)
togo.internet.simbetas <- mvrnorm(sims, togo.internet.pe, togo.internet.vc) #create simbetas
togo.internet.qoi <- logitsimfd(togo.internet.scen, togo.internet.simbetas, ci=0.95) #and get quantities of interest
togo.internet.qoi.df <- as.data.frame(togo.internet.qoi) %>% mutate(medium = 'internet', country = "togo")

#tunisia
mdata <- extractdata(internet.country, tunisia_data, na.rm = TRUE) #only extract data used in model
tunisia.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
tunisia.internet.pe <- fixef(tunisia.internet.result) # extract pe
tunisia.internet.vc <- vcov(tunisia.internet.result) # extract vcov
tunisia.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tunisia.internet.scen <- cfName(tunisia.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
tunisia.internet.scen <- cfChange(tunisia.internet.scen, "internet", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
tunisia.internet.simbetas <- mvrnorm(sims, tunisia.internet.pe, tunisia.internet.vc) #create simbetas
tunisia.internet.qoi <- logitsimfd(tunisia.internet.scen, tunisia.internet.simbetas, ci=0.95) #and get quantities of interest
tunisia.internet.qoi.df <- as.data.frame(tunisia.internet.qoi) %>% mutate(medium = 'internet', country = "tunisia")

#uganda
mdata <- extractdata(internet.country, uganda_data, na.rm = TRUE) #only extract data used in model
uganda.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
uganda.internet.pe <- fixef(uganda.internet.result) # extract pe
uganda.internet.vc <- vcov(uganda.internet.result) # extract vcov
uganda.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
uganda.internet.scen <- cfName(uganda.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
uganda.internet.scen <- cfChange(uganda.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
uganda.internet.simbetas <- mvrnorm(sims, uganda.internet.pe, uganda.internet.vc) #create simbetas
uganda.internet.qoi <- logitsimfd(uganda.internet.scen, uganda.internet.simbetas, ci=0.95) #and get quantities of interest
uganda.internet.qoi.df <- as.data.frame(uganda.internet.qoi) %>% mutate(medium = 'internet', country = "uganda")

#zambia
mdata <- extractdata(internet.country, zambia_data, na.rm = TRUE) #only extract data used in model
zambia.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
zambia.internet.pe <- fixef(zambia.internet.result) # extract pe
zambia.internet.vc <- vcov(zambia.internet.result) # extract vcov
zambia.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zambia.internet.scen <- cfName(zambia.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
zambia.internet.scen <- cfChange(zambia.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
zambia.internet.simbetas <- mvrnorm(sims, zambia.internet.pe, zambia.internet.vc) #create simbetas
zambia.internet.qoi <- logitsimfd(zambia.internet.scen, zambia.internet.simbetas, ci=0.95) #and get quantities of interest
zambia.internet.qoi.df <- as.data.frame(zambia.internet.qoi) %>% mutate(medium = 'internet', country = "zambia")

#zimbab
mdata <- extractdata(internet.country, zimbab_data, na.rm = TRUE) #only extract data used in model
zimbab.internet.result <- glmer(internet.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
zimbab.internet.pe <- fixef(zimbab.internet.result) # extract pe
zimbab.internet.vc <- vcov(zimbab.internet.result) # extract vcov
zimbab.internet.scen <- cfMake(internet.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zimbab.internet.scen <- cfName(zimbab.internet.scen, "internet 1 vs 5", scen=1) #configure scenario of 1 to 5
zimbab.internet.scen <- cfChange(zimbab.internet.scen, "internet", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
zimbab.internet.simbetas <- mvrnorm(sims, zimbab.internet.pe, zimbab.internet.vc) #create simbetas
zimbab.internet.qoi <- logitsimfd(zimbab.internet.scen, zimbab.internet.simbetas, ci=0.95) #and get quantities of interest
zimbab.internet.qoi.df <- as.data.frame(zimbab.internet.qoi) %>% mutate(medium = 'internet', country = "zimbabwe")


###############
# Social Media
###############

social_media.country <- # Use this model to estimate the result
  sexuality2 ~
  social_media +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban +
  (1 | district)  #add VSVI for country & VI for District

social_media.country.nodist <- # Use this model to create scenarios (dropping district)
  sexuality2 ~
  social_media +
  media_nosocialmedia + #consumption of other mediums
  tolerance_noLGBT2 + #add social tolerance 
  female + as.numeric(education) + as.numeric(religiosity) + as.numeric(age) + as.numeric(income_water) + urban 

#Benin
mdata <- extractdata(social_media.country, benin_data, na.rm = TRUE) #only extract data used in model
benin.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
benin.social_media.pe <- fixef(benin.social_media.result) # extract pe
benin.social_media.vc <- vcov(benin.social_media.result) # extract vcov
benin.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
benin.social_media.scen <- cfName(benin.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
benin.social_media.scen <- cfChange(benin.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
benin.social_media.simbetas <- mvrnorm(sims, benin.social_media.pe, benin.social_media.vc) #create simbetas
benin.social_media.qoi <- logitsimfd(benin.social_media.scen, benin.social_media.simbetas, ci=0.95) #and get quantities of interest
benin.social_media.qoi.df <- as.data.frame(benin.social_media.qoi) %>% mutate(medium = 'social_media', country = "benin") #create df from qoi

#Botswana
mdata <- extractdata(social_media.country.nodist, botswana_data, na.rm = TRUE) #use nodist model for botswana b/c there are no districts (i.e. no random effects)
botswana.social_media.result <- glm(social_media.country.nodist, family=binomial, data=mdata) #estimate the model
botswana.social_media.pe <-botswana.social_media.result$coefficients  # extract pe
botswana.social_media.vc <- vcov(botswana.social_media.result) # extract vcov
botswana.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
botswana.social_media.scen <- cfName(botswana.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
botswana.social_media.scen <- cfChange(botswana.social_media.scen, "social_media", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
botswana.social_media.simbetas <- mvrnorm(sims, botswana.social_media.pe, botswana.social_media.vc) #create simbetas
botswana.social_media.qoi <- logitsimfd(botswana.social_media.scen, botswana.social_media.simbetas, ci=0.95) #and get quantities of interest
botswana.social_media.qoi.df <- as.data.frame(botswana.social_media.qoi) %>% mutate(medium = 'social_media', country = "botswana")

#burkina faso
mdata <- extractdata(social_media.country, burkina_data, na.rm = TRUE) #only extract data used in model
burkina.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
burkina.social_media.pe <- fixef(burkina.social_media.result) # extract pe
burkina.social_media.vc <- vcov(burkina.social_media.result) # extract vcov
burkina.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burkina.social_media.scen <- cfName(burkina.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
burkina.social_media.scen <- cfChange(burkina.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
burkina.social_media.simbetas <- mvrnorm(sims, burkina.social_media.pe, burkina.social_media.vc) #create simbetas
burkina.social_media.qoi <- logitsimfd(burkina.social_media.scen, burkina.social_media.simbetas, ci=0.95) #and get quantities of interest
burkina.social_media.qoi.df <- as.data.frame(burkina.social_media.qoi) %>% mutate(medium = 'social_media', country = "burkina faso")

#burundi
mdata <- extractdata(social_media.country, burundi_data, na.rm = TRUE) #only extract data used in model
burundi.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
burundi.social_media.pe <- fixef(burundi.social_media.result) # extract pe
burundi.social_media.vc <- vcov(burundi.social_media.result) # extract vcov
burundi.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
burundi.social_media.scen <- cfName(burundi.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
burundi.social_media.scen <- cfChange(burundi.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
burundi.social_media.simbetas <- mvrnorm(sims, burundi.social_media.pe, burundi.social_media.vc) #create simbetas
burundi.social_media.qoi <- logitsimfd(burundi.social_media.scen, burundi.social_media.simbetas, ci=0.95) #and get quantities of interest
burundi.social_media.qoi.df <- as.data.frame(burundi.social_media.qoi) %>% mutate(medium = 'social_media', country = "burundi")

#cameroon
mdata <- extractdata(social_media.country, cameroon_data, na.rm = TRUE) #only extract data used in model
cameroon.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                  nAGQ = 0 ) # optimizer for faster but less precise fit
cameroon.social_media.pe <- fixef(cameroon.social_media.result) # extract pe
cameroon.social_media.vc <- vcov(cameroon.social_media.result) # extract vcov
cameroon.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cameroon.social_media.scen <- cfName(cameroon.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
cameroon.social_media.scen <- cfChange(cameroon.social_media.scen, "social_media", 
                                   x = 5,
                                   xpre=1,
                                   scen=1)
cameroon.social_media.simbetas <- mvrnorm(sims, cameroon.social_media.pe, cameroon.social_media.vc) #create simbetas
cameroon.social_media.qoi <- logitsimfd(cameroon.social_media.scen, cameroon.social_media.simbetas, ci=0.95) #and get quantities of interest
cameroon.social_media.qoi.df <- as.data.frame(cameroon.social_media.qoi) %>% mutate(medium = 'social_media', country = "cameroon")

#cape verde
mdata <- extractdata(social_media.country, capeverde_data, na.rm = TRUE) #only extract data used in model
capeverde.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                   nAGQ = 0 ) # optimizer for faster but less precise fit
capeverde.social_media.pe <- fixef(capeverde.social_media.result) # extract pe
capeverde.social_media.vc <- vcov(capeverde.social_media.result) # extract vcov
capeverde.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
capeverde.social_media.scen <- cfName(capeverde.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
capeverde.social_media.scen <- cfChange(capeverde.social_media.scen, "social_media", 
                                    x = 5,
                                    xpre=1,
                                    scen=1)
capeverde.social_media.simbetas <- mvrnorm(sims, capeverde.social_media.pe, capeverde.social_media.vc) #create simbetas
capeverde.social_media.qoi <- logitsimfd(capeverde.social_media.scen, capeverde.social_media.simbetas, ci=0.95) #and get quantities of interest
capeverde.social_media.qoi.df <- as.data.frame(capeverde.social_media.qoi) %>% mutate(medium = 'social_media', country = "cape verde")

#cote d'ivoire
mdata <- extractdata(social_media.country, cdi_data, na.rm = TRUE) #only extract data used in model
cdi.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                             nAGQ = 0 ) # optimizer for faster but less precise fit
cdi.social_media.pe <- fixef(cdi.social_media.result) # extract pe
cdi.social_media.vc <- vcov(cdi.social_media.result) # extract vcov
cdi.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
cdi.social_media.scen <- cfName(cdi.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
cdi.social_media.scen <- cfChange(cdi.social_media.scen, "social_media", 
                              x = 5,
                              xpre=1,
                              scen=1)
cdi.social_media.simbetas <- mvrnorm(sims, cdi.social_media.pe, cdi.social_media.vc) #create simbetas
cdi.social_media.qoi <- logitsimfd(cdi.social_media.scen, cdi.social_media.simbetas, ci=0.95) #and get quantities of interest
cdi.social_media.qoi.df <- as.data.frame(cdi.social_media.qoi) %>% mutate(medium = 'social_media', country = "cote d'ivoire")

#gabon
mdata <- extractdata(social_media.country, gabon_data, na.rm = TRUE) #only extract data used in model
gabon.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
gabon.social_media.pe <- fixef(gabon.social_media.result) # extract pe
gabon.social_media.vc <- vcov(gabon.social_media.result) # extract vcov
gabon.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
gabon.social_media.scen <- cfName(gabon.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
gabon.social_media.scen <- cfChange(gabon.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
gabon.social_media.simbetas <- mvrnorm(sims, gabon.social_media.pe, gabon.social_media.vc) #create simbetas
gabon.social_media.qoi <- logitsimfd(gabon.social_media.scen, gabon.social_media.simbetas, ci=0.95) #and get quantities of interest
gabon.social_media.qoi.df <- as.data.frame(gabon.social_media.qoi) %>% mutate(medium = 'social_media', country = "gabon")

#ghana
mdata <- extractdata(social_media.country, ghana_data, na.rm = TRUE) #only extract data used in model
ghana.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
ghana.social_media.pe <- fixef(ghana.social_media.result) # extract pe
ghana.social_media.vc <- vcov(ghana.social_media.result) # extract vcov
ghana.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
ghana.social_media.scen <- cfName(ghana.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
ghana.social_media.scen <- cfChange(ghana.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
ghana.social_media.simbetas <- mvrnorm(sims, ghana.social_media.pe, ghana.social_media.vc) #create simbetas
ghana.social_media.qoi <- logitsimfd(ghana.social_media.scen, ghana.social_media.simbetas, ci=0.95) #and get quantities of interest
ghana.social_media.qoi.df <- as.data.frame(ghana.social_media.qoi) %>% mutate(medium = 'social_media', country = "ghana")

#guinea 
mdata <- extractdata(social_media.country, guinea_data, na.rm = TRUE) #only extract data used in model
guinea.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
guinea.social_media.pe <- fixef(guinea.social_media.result) # extract pe
guinea.social_media.vc <- vcov(guinea.social_media.result) # extract vcov
guinea.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
guinea.social_media.scen <- cfName(guinea.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
guinea.social_media.scen <- cfChange(guinea.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
guinea.social_media.simbetas <- mvrnorm(sims, guinea.social_media.pe, guinea.social_media.vc) #create simbetas
guinea.social_media.qoi <- logitsimfd(guinea.social_media.scen, guinea.social_media.simbetas, ci=0.95) #and get quantities of interest
guinea.social_media.qoi.df <- as.data.frame(guinea.social_media.qoi) %>% mutate(medium = 'social_media', country = "guinea")

#kenya
mdata <- extractdata(social_media.country, kenya_data, na.rm = TRUE) #only extract data used in model
kenya.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
kenya.social_media.pe <- fixef(kenya.social_media.result) # extract pe
kenya.social_media.vc <- vcov(kenya.social_media.result) # extract vcov
kenya.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
kenya.social_media.scen <- cfName(kenya.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
kenya.social_media.scen <- cfChange(kenya.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
kenya.social_media.simbetas <- mvrnorm(sims, kenya.social_media.pe, kenya.social_media.vc) #create simbetas
kenya.social_media.qoi <- logitsimfd(kenya.social_media.scen, kenya.social_media.simbetas, ci=0.95) #and get quantities of interest
kenya.social_media.qoi.df <- as.data.frame(kenya.social_media.qoi) %>% mutate(medium = 'social_media', country = "kenya")

#lesotho
mdata <- extractdata(social_media.country, lesotho_data, na.rm = TRUE) #only extract data used in model
lesotho.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
lesotho.social_media.pe <- fixef(lesotho.social_media.result) # extract pe
lesotho.social_media.vc <- vcov(lesotho.social_media.result) # extract vcov
lesotho.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
lesotho.social_media.scen <- cfName(lesotho.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
lesotho.social_media.scen <- cfChange(lesotho.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
lesotho.social_media.simbetas <- mvrnorm(sims, lesotho.social_media.pe, lesotho.social_media.vc) #create simbetas
lesotho.social_media.qoi <- logitsimfd(lesotho.social_media.scen, lesotho.social_media.simbetas, ci=0.95) #and get quantities of interest
lesotho.social_media.qoi.df <- as.data.frame(lesotho.social_media.qoi) %>% mutate(medium = 'social_media', country = "lesotho")

#liberia
mdata <- extractdata(social_media.country, liberia_data, na.rm = TRUE) #only extract data used in model
liberia.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
liberia.social_media.pe <- fixef(liberia.social_media.result) # extract pe
liberia.social_media.vc <- vcov(liberia.social_media.result) # extract vcov
liberia.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
liberia.social_media.scen <- cfName(liberia.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
liberia.social_media.scen <- cfChange(liberia.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
liberia.social_media.simbetas <- mvrnorm(sims, liberia.social_media.pe, liberia.social_media.vc) #create simbetas
liberia.social_media.qoi <- logitsimfd(liberia.social_media.scen, liberia.social_media.simbetas, ci=0.95) #and get quantities of interest
liberia.social_media.qoi.df <- as.data.frame(liberia.social_media.qoi) %>% mutate(medium = 'social_media', country = "liberia")

#madagascar
mdata <- extractdata(social_media.country, madag_data, na.rm = TRUE) #only extract data used in model
madag.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
madag.social_media.pe <- fixef(madag.social_media.result) # extract pe
madag.social_media.vc <- vcov(madag.social_media.result) # extract vcov
madag.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
madag.social_media.scen <- cfName(madag.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
madag.social_media.scen <- cfChange(madag.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
madag.social_media.simbetas <- mvrnorm(sims, madag.social_media.pe, madag.social_media.vc) #create simbetas
madag.social_media.qoi <- logitsimfd(madag.social_media.scen, madag.social_media.simbetas, ci=0.95) #and get quantities of interest
madag.social_media.qoi.df <- as.data.frame(madag.social_media.qoi) %>% mutate(medium = 'social_media', country = "madagascar")

#malawi
mdata <- extractdata(social_media.country, malawi_data, na.rm = TRUE) #only extract data used in model
malawi.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
malawi.social_media.pe <- fixef(malawi.social_media.result) # extract pe
malawi.social_media.vc <- vcov(malawi.social_media.result) # extract vcov
malawi.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
malawi.social_media.scen <- cfName(malawi.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
malawi.social_media.scen <- cfChange(malawi.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
malawi.social_media.simbetas <- mvrnorm(sims, malawi.social_media.pe, malawi.social_media.vc) #create simbetas
malawi.social_media.qoi <- logitsimfd(malawi.social_media.scen, malawi.social_media.simbetas, ci=0.95) #and get quantities of interest
malawi.social_media.qoi.df <- as.data.frame(malawi.social_media.qoi) %>% mutate(medium = 'social_media', country = "malawi")

#mali
mdata <- extractdata(social_media.country, mali_data, na.rm = TRUE) #only extract data used in model
mali.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
mali.social_media.pe <- fixef(mali.social_media.result) # extract pe
mali.social_media.vc <- vcov(mali.social_media.result) # extract vcov
mali.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mali.social_media.scen <- cfName(mali.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
mali.social_media.scen <- cfChange(mali.social_media.scen, "social_media", 
                               x = 5,
                               xpre=1,
                               scen=1)
mali.social_media.simbetas <- mvrnorm(sims, mali.social_media.pe, mali.social_media.vc) #create simbetas
mali.social_media.qoi <- logitsimfd(mali.social_media.scen, mali.social_media.simbetas, ci=0.95) #and get quantities of interest
mali.social_media.qoi.df <- as.data.frame(mali.social_media.qoi) %>% mutate(medium = 'social_media', country = "mali")

#mauritius
mdata <- extractdata(social_media.country, maurit_data, na.rm = TRUE) #only extract data used in model
maurit.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
maurit.social_media.pe <- fixef(maurit.social_media.result) # extract pe
maurit.social_media.vc <- vcov(maurit.social_media.result) # extract vcov
maurit.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
maurit.social_media.scen <- cfName(maurit.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
maurit.social_media.scen <- cfChange(maurit.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
maurit.social_media.simbetas <- mvrnorm(sims, maurit.social_media.pe, maurit.social_media.vc) #create simbetas
maurit.social_media.qoi <- logitsimfd(maurit.social_media.scen, maurit.social_media.simbetas, ci=0.95) #and get quantities of interest
maurit.social_media.qoi.df <- as.data.frame(maurit.social_media.qoi) %>% mutate(medium = 'social_media', country = "mauritius")

#morocco
mdata <- extractdata(social_media.country, morocco_data, na.rm = TRUE) #only extract data used in model
morocco.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
morocco.social_media.pe <- fixef(morocco.social_media.result) # extract pe
morocco.social_media.vc <- vcov(morocco.social_media.result) # extract vcov
morocco.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
morocco.social_media.scen <- cfName(morocco.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
morocco.social_media.scen <- cfChange(morocco.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
morocco.social_media.simbetas <- mvrnorm(sims, morocco.social_media.pe, morocco.social_media.vc) #create simbetas
morocco.social_media.qoi <- logitsimfd(morocco.social_media.scen, morocco.social_media.simbetas, ci=0.95) #and get quantities of interest
morocco.social_media.qoi.df <- as.data.frame(morocco.social_media.qoi) %>% mutate(medium = 'social_media', country = "morocco")

#mozambique
mdata <- extractdata(social_media.country, mozamb_data, na.rm = TRUE) #only extract data used in model
mozamb.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
mozamb.social_media.pe <- fixef(mozamb.social_media.result) # extract pe
mozamb.social_media.vc <- vcov(mozamb.social_media.result) # extract vcov
mozamb.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
mozamb.social_media.scen <- cfName(mozamb.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
mozamb.social_media.scen <- cfChange(mozamb.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
mozamb.social_media.simbetas <- mvrnorm(sims, mozamb.social_media.pe, mozamb.social_media.vc) #create simbetas
mozamb.social_media.qoi <- logitsimfd(mozamb.social_media.scen, mozamb.social_media.simbetas, ci=0.95) #and get quantities of interest
mozamb.social_media.qoi.df <- as.data.frame(mozamb.social_media.qoi) %>% mutate(medium = 'social_media', country = "mozambique")

#namibia
mdata <- extractdata(social_media.country, namibia_data, na.rm = TRUE) #only extract data used in model
namibia.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
namibia.social_media.pe <- fixef(namibia.social_media.result) # extract pe
namibia.social_media.vc <- vcov(namibia.social_media.result) # extract vcov
namibia.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
namibia.social_media.scen <- cfName(namibia.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
namibia.social_media.scen <- cfChange(namibia.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
namibia.social_media.simbetas <- mvrnorm(sims, namibia.social_media.pe, namibia.social_media.vc) #create simbetas
namibia.social_media.qoi <- logitsimfd(namibia.social_media.scen, namibia.social_media.simbetas, ci=0.95) #and get quantities of interest
namibia.social_media.qoi.df <- as.data.frame(namibia.social_media.qoi) %>% mutate(medium = 'social_media', country = "namibia")

#niger
mdata <- extractdata(social_media.country, niger_data, na.rm = TRUE) #only extract data used in model
niger.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                               nAGQ = 0 ) # optimizer for faster but less precise fit
niger.social_media.pe <- fixef(niger.social_media.result) # extract pe
niger.social_media.vc <- vcov(niger.social_media.result) # extract vcov
niger.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
niger.social_media.scen <- cfName(niger.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
niger.social_media.scen <- cfChange(niger.social_media.scen, "social_media", 
                                x = 5,
                                xpre=1,
                                scen=1)
niger.social_media.simbetas <- mvrnorm(sims, niger.social_media.pe, niger.social_media.vc) #create simbetas
niger.social_media.qoi <- logitsimfd(niger.social_media.scen, niger.social_media.simbetas, ci=0.95) #and get quantities of interest
niger.social_media.qoi.df <- as.data.frame(niger.social_media.qoi) %>% mutate(medium = 'social_media', country = "niger")

#nigeria
mdata <- extractdata(social_media.country, nigeria_data, na.rm = TRUE) #only extract data used in model
nigeria.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
nigeria.social_media.pe <- fixef(nigeria.social_media.result) # extract pe
nigeria.social_media.vc <- vcov(nigeria.social_media.result) # extract vcov
nigeria.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
nigeria.social_media.scen <- cfName(nigeria.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
nigeria.social_media.scen <- cfChange(nigeria.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
nigeria.social_media.simbetas <- mvrnorm(sims, nigeria.social_media.pe, nigeria.social_media.vc) #create simbetas
nigeria.social_media.qoi <- logitsimfd(nigeria.social_media.scen, nigeria.social_media.simbetas, ci=0.95) #and get quantities of interest
nigeria.social_media.qoi.df <- as.data.frame(nigeria.social_media.qoi) %>% mutate(medium = 'social_media', country = "nigeria")

#são tomé and príncipe
mdata <- extractdata(social_media.country, stp_data, na.rm = TRUE) #only extract data used in model
stp.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                             nAGQ = 0 ) # optimizer for faster but less precise fit
stp.social_media.pe <- fixef(stp.social_media.result) # extract pe
stp.social_media.vc <- vcov(stp.social_media.result) # extract vcov
stp.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
stp.social_media.scen <- cfName(stp.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
stp.social_media.scen <- cfChange(stp.social_media.scen, "social_media", 
                              x = 5,
                              xpre=1,
                              scen=1)
stp.social_media.simbetas <- mvrnorm(sims, stp.social_media.pe, stp.social_media.vc) #create simbetas
stp.social_media.qoi <- logitsimfd(stp.social_media.scen, stp.social_media.simbetas, ci=0.95) #and get quantities of interest
stp.social_media.qoi.df <- as.data.frame(stp.social_media.qoi) %>% mutate(medium = 'social_media', country = "são tomé and príncipe")

#senegal
mdata <- extractdata(social_media.country, senegal_data, na.rm = TRUE) #only extract data used in model
senegal.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
senegal.social_media.pe <- fixef(senegal.social_media.result) # extract pe
senegal.social_media.vc <- vcov(senegal.social_media.result) # extract vcov
senegal.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
senegal.social_media.scen <- cfName(senegal.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
senegal.social_media.scen <- cfChange(senegal.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
senegal.social_media.simbetas <- mvrnorm(sims, senegal.social_media.pe, senegal.social_media.vc) #create simbetas
senegal.social_media.qoi <- logitsimfd(senegal.social_media.scen, senegal.social_media.simbetas, ci=0.95) #and get quantities of interest
senegal.social_media.qoi.df <- as.data.frame(senegal.social_media.qoi) %>% mutate(medium = 'social_media', country = "senegal")

#sierra
mdata <- extractdata(social_media.country, sierra_data, na.rm = TRUE) #only extract data used in model
sierra.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
sierra.social_media.pe <- fixef(sierra.social_media.result) # extract pe
sierra.social_media.vc <- vcov(sierra.social_media.result) # extract vcov
sierra.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
sierra.social_media.scen <- cfName(sierra.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
sierra.social_media.scen <- cfChange(sierra.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
sierra.social_media.simbetas <- mvrnorm(sims, sierra.social_media.pe, sierra.social_media.vc) #create simbetas
sierra.social_media.qoi <- logitsimfd(sierra.social_media.scen, sierra.social_media.simbetas, ci=0.95) #and get quantities of interest
sierra.social_media.qoi.df <- as.data.frame(sierra.social_media.qoi) %>% mutate(medium = 'social_media', country = "sierra leone")

#southaf
mdata <- extractdata(social_media.country, southaf_data, na.rm = TRUE) #only extract data used in model
southaf.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
southaf.social_media.pe <- fixef(southaf.social_media.result) # extract pe
southaf.social_media.vc <- vcov(southaf.social_media.result) # extract vcov
southaf.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
southaf.social_media.scen <- cfName(southaf.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
southaf.social_media.scen <- cfChange(southaf.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
southaf.social_media.simbetas <- mvrnorm(sims, southaf.social_media.pe, southaf.social_media.vc) #create simbetas
southaf.social_media.qoi <- logitsimfd(southaf.social_media.scen, southaf.social_media.simbetas, ci=0.95) #and get quantities of interest
southaf.social_media.qoi.df <- as.data.frame(southaf.social_media.qoi) %>% mutate(medium = 'social_media', country = "south africa")

#swaz
mdata <- extractdata(social_media.country, swaz_data, na.rm = TRUE) #only extract data used in model
swaz.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
swaz.social_media.pe <- fixef(swaz.social_media.result) # extract pe
swaz.social_media.vc <- vcov(swaz.social_media.result) # extract vcov
swaz.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
swaz.social_media.scen <- cfName(swaz.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
swaz.social_media.scen <- cfChange(swaz.social_media.scen, "social_media", 
                               x = 5,
                               xpre=1,
                               scen=1)
swaz.social_media.simbetas <- mvrnorm(sims, swaz.social_media.pe, swaz.social_media.vc) #create simbetas
swaz.social_media.qoi <- logitsimfd(swaz.social_media.scen, swaz.social_media.simbetas, ci=0.95) #and get quantities of interest
swaz.social_media.qoi.df <- as.data.frame(swaz.social_media.qoi) %>% mutate(medium = 'social_media', country = "swaziland")

#tanz
mdata <- extractdata(social_media.country, tanz_data, na.rm = TRUE) #only extract data used in model
tanz.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
tanz.social_media.pe <- fixef(tanz.social_media.result) # extract pe
tanz.social_media.vc <- vcov(tanz.social_media.result) # extract vcov
tanz.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tanz.social_media.scen <- cfName(tanz.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
tanz.social_media.scen <- cfChange(tanz.social_media.scen, "social_media", 
                               x = 5,
                               xpre=1,
                               scen=1)
tanz.social_media.simbetas <- mvrnorm(sims, tanz.social_media.pe, tanz.social_media.vc) #create simbetas
tanz.social_media.qoi <- logitsimfd(tanz.social_media.scen, tanz.social_media.simbetas, ci=0.95) #and get quantities of interest
tanz.social_media.qoi.df <- as.data.frame(tanz.social_media.qoi) %>% mutate(medium = 'social_media', country = "tanzania")

#togo
mdata <- extractdata(social_media.country, togo_data, na.rm = TRUE) #only extract data used in model
togo.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                              nAGQ = 0 ) # optimizer for faster but less precise fit
togo.social_media.pe <- fixef(togo.social_media.result) # extract pe
togo.social_media.vc <- vcov(togo.social_media.result) # extract vcov
togo.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
togo.social_media.scen <- cfName(togo.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
togo.social_media.scen <- cfChange(togo.social_media.scen, "social_media", 
                               x = 5,
                               xpre=1,
                               scen=1)
togo.social_media.simbetas <- mvrnorm(sims, togo.social_media.pe, togo.social_media.vc) #create simbetas
togo.social_media.qoi <- logitsimfd(togo.social_media.scen, togo.social_media.simbetas, ci=0.95) #and get quantities of interest
togo.social_media.qoi.df <- as.data.frame(togo.social_media.qoi) %>% mutate(medium = 'social_media', country = "togo")

#tunisia
mdata <- extractdata(social_media.country, tunisia_data, na.rm = TRUE) #only extract data used in model
tunisia.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                 nAGQ = 0 ) # optimizer for faster but less precise fit
tunisia.social_media.pe <- fixef(tunisia.social_media.result) # extract pe
tunisia.social_media.vc <- vcov(tunisia.social_media.result) # extract vcov
tunisia.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
tunisia.social_media.scen <- cfName(tunisia.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
tunisia.social_media.scen <- cfChange(tunisia.social_media.scen, "social_media", 
                                  x = 5,
                                  xpre=1,
                                  scen=1)
tunisia.social_media.simbetas <- mvrnorm(sims, tunisia.social_media.pe, tunisia.social_media.vc) #create simbetas
tunisia.social_media.qoi <- logitsimfd(tunisia.social_media.scen, tunisia.social_media.simbetas, ci=0.95) #and get quantities of interest
tunisia.social_media.qoi.df <- as.data.frame(tunisia.social_media.qoi) %>% mutate(medium = 'social_media', country = "tunisia")

#uganda
mdata <- extractdata(social_media.country, uganda_data, na.rm = TRUE) #only extract data used in model
uganda.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
uganda.social_media.pe <- fixef(uganda.social_media.result) # extract pe
uganda.social_media.vc <- vcov(uganda.social_media.result) # extract vcov
uganda.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
uganda.social_media.scen <- cfName(uganda.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
uganda.social_media.scen <- cfChange(uganda.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
uganda.social_media.simbetas <- mvrnorm(sims, uganda.social_media.pe, uganda.social_media.vc) #create simbetas
uganda.social_media.qoi <- logitsimfd(uganda.social_media.scen, uganda.social_media.simbetas, ci=0.95) #and get quantities of interest
uganda.social_media.qoi.df <- as.data.frame(uganda.social_media.qoi) %>% mutate(medium = 'social_media', country = "uganda")

#zambia
mdata <- extractdata(social_media.country, zambia_data, na.rm = TRUE) #only extract data used in model
zambia.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
zambia.social_media.pe <- fixef(zambia.social_media.result) # extract pe
zambia.social_media.vc <- vcov(zambia.social_media.result) # extract vcov
zambia.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zambia.social_media.scen <- cfName(zambia.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
zambia.social_media.scen <- cfChange(zambia.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
zambia.social_media.simbetas <- mvrnorm(sims, zambia.social_media.pe, zambia.social_media.vc) #create simbetas
zambia.social_media.qoi <- logitsimfd(zambia.social_media.scen, zambia.social_media.simbetas, ci=0.95) #and get quantities of interest
zambia.social_media.qoi.df <- as.data.frame(zambia.social_media.qoi) %>% mutate(medium = 'social_media', country = "zambia")

#zimbab
mdata <- extractdata(social_media.country, zimbab_data, na.rm = TRUE) #only extract data used in model
zimbab.social_media.result <- glmer(social_media.country, family=binomial, data=mdata, #estimate the model
                                nAGQ = 0 ) # optimizer for faster but less precise fit
zimbab.social_media.pe <- fixef(zimbab.social_media.result) # extract pe
zimbab.social_media.vc <- vcov(zimbab.social_media.result) # extract vcov
zimbab.social_media.scen <- cfMake(social_media.country.nodist, data=mdata, nscen=1) #initialize scenario to mean value of covariates
zimbab.social_media.scen <- cfName(zimbab.social_media.scen, "social_media 1 vs 5", scen=1) #configure scenario of 1 to 5
zimbab.social_media.scen <- cfChange(zimbab.social_media.scen, "social_media", 
                                 x = 5,
                                 xpre=1,
                                 scen=1)
zimbab.social_media.simbetas <- mvrnorm(sims, zimbab.social_media.pe, zimbab.social_media.vc) #create simbetas
zimbab.social_media.qoi <- logitsimfd(zimbab.social_media.scen, zimbab.social_media.simbetas, ci=0.95) #and get quantities of interest
zimbab.social_media.qoi.df <- as.data.frame(zimbab.social_media.qoi) %>% mutate(medium = 'social_media', country = "zimbabwe")

########################
# Combine in dataframe
# & plot
########################

# a dataframe of FD RADIO results for every country
radio.fd.df <- as.data.frame(bind_rows(benin.radio.qoi.df, botswana.radio.qoi.df, burkina.radio.qoi.df, burundi.radio.qoi.df, cameroon.radio.qoi.df,
                                       capeverde.radio.qoi.df, cdi.radio.qoi.df, gabon.radio.qoi.df, ghana.radio.qoi.df, guinea.radio.qoi.df,
                                       kenya.radio.qoi.df, lesotho.radio.qoi.df, liberia.radio.qoi.df, madag.radio.qoi.df, malawi.radio.qoi.df,
                                       mali.radio.qoi.df, maurit.radio.qoi.df, morocco.radio.qoi.df, mozamb.radio.qoi.df, namibia.radio.qoi.df, 
                                       niger.radio.qoi.df, nigeria.radio.qoi.df, stp.radio.qoi.df, senegal.radio.qoi.df, sierra.radio.qoi.df,
                                       southaf.radio.qoi.df, swaz.radio.qoi.df, tanz.radio.qoi.df, togo.radio.qoi.df, tunisia.radio.qoi.df, 
                                       uganda.radio.qoi.df, zambia.radio.qoi.df, zimbab.radio.qoi.df))
# a dataframe of FD TV results for every country
tv.fd.df <- as.data.frame(bind_rows(benin.tv.qoi.df, botswana.tv.qoi.df, burkina.tv.qoi.df, burundi.tv.qoi.df, cameroon.tv.qoi.df,
                                       capeverde.tv.qoi.df, cdi.tv.qoi.df, gabon.tv.qoi.df, ghana.tv.qoi.df, guinea.tv.qoi.df,
                                       kenya.tv.qoi.df, lesotho.tv.qoi.df, liberia.tv.qoi.df, madag.tv.qoi.df, malawi.tv.qoi.df,
                                       mali.tv.qoi.df, maurit.tv.qoi.df, morocco.tv.qoi.df, mozamb.tv.qoi.df, namibia.tv.qoi.df, 
                                       niger.tv.qoi.df, nigeria.tv.qoi.df, stp.tv.qoi.df, senegal.tv.qoi.df, sierra.tv.qoi.df,
                                       southaf.tv.qoi.df, swaz.tv.qoi.df, tanz.tv.qoi.df, togo.tv.qoi.df, tunisia.tv.qoi.df, 
                                       uganda.tv.qoi.df, zambia.tv.qoi.df, zimbab.tv.qoi.df))
# a dataframe of FD NEWSPAPER results for every country
newspaper.fd.df <- as.data.frame(bind_rows(benin.newspaper.qoi.df, botswana.newspaper.qoi.df, burkina.newspaper.qoi.df, burundi.newspaper.qoi.df, cameroon.newspaper.qoi.df,
                                    capeverde.newspaper.qoi.df, cdi.newspaper.qoi.df, gabon.newspaper.qoi.df, ghana.newspaper.qoi.df, guinea.newspaper.qoi.df,
                                    kenya.newspaper.qoi.df, lesotho.newspaper.qoi.df, liberia.newspaper.qoi.df, madag.newspaper.qoi.df, malawi.newspaper.qoi.df,
                                    mali.newspaper.qoi.df, maurit.newspaper.qoi.df, morocco.newspaper.qoi.df, mozamb.newspaper.qoi.df, namibia.newspaper.qoi.df, 
                                    niger.newspaper.qoi.df, nigeria.newspaper.qoi.df, stp.newspaper.qoi.df, senegal.newspaper.qoi.df, sierra.newspaper.qoi.df,
                                    southaf.newspaper.qoi.df, swaz.newspaper.qoi.df, tanz.newspaper.qoi.df, togo.newspaper.qoi.df, tunisia.newspaper.qoi.df, 
                                    uganda.newspaper.qoi.df, zambia.newspaper.qoi.df, zimbab.newspaper.qoi.df))
# a dataframe of FD INTERNET results for every country
internet.fd.df <- as.data.frame(bind_rows(benin.internet.qoi.df, botswana.internet.qoi.df, burkina.internet.qoi.df, burundi.internet.qoi.df, cameroon.internet.qoi.df,
                                    capeverde.internet.qoi.df, cdi.internet.qoi.df, gabon.internet.qoi.df, ghana.internet.qoi.df, guinea.internet.qoi.df,
                                    kenya.internet.qoi.df, lesotho.internet.qoi.df, liberia.internet.qoi.df, madag.internet.qoi.df, malawi.internet.qoi.df,
                                    mali.internet.qoi.df, maurit.internet.qoi.df, morocco.internet.qoi.df, mozamb.internet.qoi.df, namibia.internet.qoi.df, 
                                    niger.internet.qoi.df, nigeria.internet.qoi.df, stp.internet.qoi.df, senegal.internet.qoi.df, sierra.internet.qoi.df,
                                    southaf.internet.qoi.df, swaz.internet.qoi.df, tanz.internet.qoi.df, togo.internet.qoi.df, tunisia.internet.qoi.df, 
                                    uganda.internet.qoi.df, zambia.internet.qoi.df, zimbab.internet.qoi.df))
# a dataframe of FD SOCIAL MEDIA results for every country
social_media.fd.df <- as.data.frame(bind_rows(benin.social_media.qoi.df, botswana.social_media.qoi.df, burkina.social_media.qoi.df, burundi.social_media.qoi.df, cameroon.social_media.qoi.df,
                                    capeverde.social_media.qoi.df, cdi.social_media.qoi.df, gabon.social_media.qoi.df, ghana.social_media.qoi.df, guinea.social_media.qoi.df,
                                    kenya.social_media.qoi.df, lesotho.social_media.qoi.df, liberia.social_media.qoi.df, madag.social_media.qoi.df, malawi.social_media.qoi.df,
                                    mali.social_media.qoi.df, maurit.social_media.qoi.df, morocco.social_media.qoi.df, mozamb.social_media.qoi.df, namibia.social_media.qoi.df, 
                                    niger.social_media.qoi.df, nigeria.social_media.qoi.df, stp.social_media.qoi.df, senegal.social_media.qoi.df, sierra.social_media.qoi.df,
                                    southaf.social_media.qoi.df, swaz.social_media.qoi.df, tanz.social_media.qoi.df, togo.social_media.qoi.df, tunisia.social_media.qoi.df, 
                                    uganda.social_media.qoi.df, zambia.social_media.qoi.df, zimbab.social_media.qoi.df))

# add KOF rating to each country
kof_africa <- read.csv("data-raw/kof/Data_2018.csv") #read raw data
kof_africa$country <- as.character(kof_africa$country) # rename country names to match Afrob
kof_africa$country[kof_africa$country == 'Egypt, Arab Rep.'] <- "Egypt"
kof_africa$country[kof_africa$country == 'Sao Tome and Principe'] <- "São Tomé and Príncipe"
afrob_names <-    # create vector of Afrob country names
  as.vector(unique(myData$country)) 
kof_africa <- kof_africa %>%    # subset to keep country in Afrob and latest year of data
  subset(country %in% afrob_names) %>%
  subset(year == "2015") %>%
  dplyr::select(country, KOFSoGI) #keep only the KOF score I'm using

kof_africa$country <- tolower(kof_africa$country)
radio.fd.df <- left_join(radio.fd.df, kof_africa, by = "country")
tv.fd.df <- left_join(tv.fd.df, kof_africa, by = "country")
newspaper.fd.df <- left_join(newspaper.fd.df, kof_africa, by = "country")
internet.fd.df <- left_join(internet.fd.df, kof_africa, by = "country")
social_media.fd.df <- left_join(social_media.fd.df, kof_africa, by = "country")

# PLOT w/ KOF SCORE
# a plot of Radio FD...
# Figure A.6a
radio_country_kof_fd_plot <- ggplot(radio.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, KOFSoGI) #reorder based on KOF score
  )) +
  labs(x = "country (listed in descending order by KOF score)", y = "difference in probability of supporting LGBTQs (95% CI)") +
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
radio_country_kof_fd_plot

ggsave(filename = "figures/radio_country_kof_ml_fd_plot.pdf", #Save it to directory
       plot = radio_country_kof_fd_plot,
       width = 8, height = 6)

# a plot of TV FD ...
# Figure A.6b
tv_country_kof_fd_plot <- ggplot(tv.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, KOFSoGI) #reorder based on KOF score
  )) +
  labs(x = "country (listed in descending order by KOF score)", y = "difference in probability of supporting LGBTQs (95% CI)") +
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
tv_country_kof_fd_plot

ggsave(filename = "figures/tv_country_kof_ml_fd_plot.pdf", #Save it to directory
       plot = tv_country_kof_fd_plot,
       width = 8, height = 6)

# a plot of NEWSPAPER FD ...
# Figure A.6c
newspaper_country_kof_fd_plot <- ggplot(newspaper.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, KOFSoGI) #reorder based on KOF score
  )) +
  labs(x = "country (listed in descending order by KOF score)", y = "difference in probability of supporting LGBTQs (95% CI)") +
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray80", "gray50")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
newspaper_country_kof_fd_plot

ggsave(filename = "figures/newspaper_country_kof_ml_fd_plot.pdf", #Save it to directory
       plot = newspaper_country_kof_fd_plot,
       width = 8, height = 6)

# a plot of INTERNET FD ...
# Figure A.6d
internet_country_kof_fd_plot <- ggplot(internet.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, KOFSoGI) #reorder based on fh score
  )) +
  labs(x = "country (listed in descending order by KOF score)", y = "difference in probability of supporting LGBTQs (95% CI)") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' internet consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
internet_country_kof_fd_plot

ggsave(filename = "figures/internet_country_kof_ml_fd_plot.pdf", #Save it to directory
       plot = internet_country_kof_fd_plot,
       width = 8, height = 6)

# a plot of SOCIAL MEDIA FD ...
# Figure A.6e
social_media_country_kof_fd_plot <- ggplot(social_media.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, KOFSoGI) #reorder based on KOF score
  )) +
  labs(x = "country (listed in descending order by KOF score)", y = "difference in probability of supporting LGBTQs (95% CI)") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' social media consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
social_media_country_kof_fd_plot

ggsave(filename = "figures/social_media_country_kof_ml_fd_plot.pdf", #Save it to directory
       plot = social_media_country_kof_fd_plot,
       width = 8, height = 6)

# add FH status to each country
fh_africa <- readRDS("data/freedom_house_africa.rds")
fh_africa$fh_scale <- as.numeric(fh_africa$fh_scale) #convert to numeric so I can order it in ggplot
fh_africa$country <- tolower(fh_africa$country)
radio.fd.df <- left_join(radio.fd.df, fh_africa, by = "country")
tv.fd.df <- left_join(tv.fd.df, fh_africa, by = "country")
newspaper.fd.df <- left_join(newspaper.fd.df, fh_africa, by = "country")
internet.fd.df <- left_join(internet.fd.df, fh_africa, by = "country")
social_media.fd.df <- left_join(social_media.fd.df, fh_africa, by = "country")

# PLOT w/ FH STATUS
# a plot of Radio FD...
# Figure A.5a
radio.fd.df %<>% subset(!is.na(fh_scale)) # remove sao tome b/c it doesn't have FH score
radio.fd.df %<>% mutate(fh_status = fct_relevel(fh_status, #relevel so the legend is in order
                                                c("F", "PF", "NF")))
radio_country_fd_plot <- ggplot(radio.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, fh_scale), #reorder based on fh score
                      color=fh_status
                      )) +
  labs(x = "country (listed in descending order by Freedom House score)", y = "difference in probability of supporting LGBTQs (95% CI)",
       caption = "F = Free, PF = Partly Free, NF = Not Free") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' radio consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray50", "gray80")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
                     )
radio_country_fd_plot

ggsave(filename = "figures/radio_country_ml_fd_plot.pdf", #Save it to directory
       plot = radio_country_fd_plot,
       width = 8, height = 6)

# a plot of TV FD ...
# Figure A.5b
tv.fd.df %<>% subset(!is.na(fh_scale)) # remove sao tome b/c it doesn't have FH score
tv.fd.df %<>% mutate(fh_status = fct_relevel(fh_status, #relevel so the legend is in order
                                                c("F", "PF", "NF")))
tv_country_fd_plot <- ggplot(tv.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, fh_scale), #reorder based on fh score
                      color=fh_status
  )) +
  labs(x = "country (listed in descending order by Freedom House score)", y = "difference in probability of supporting LGBTQs (95% CI)",
       caption = "F = Free, PF = Partly Free, NF = Not Free") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' tv consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray50", "gray80")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
tv_country_fd_plot

ggsave(filename = "figures/tv_country_ml_fd_plot.pdf", #Save it to directory
       plot = tv_country_fd_plot,
       width = 8, height = 6)


# a plot of NEWSPAPER FD ...
# Figure A.5c
newspaper.fd.df %<>% subset(!is.na(fh_scale)) # remove sao tome b/c it doesn't have FH score
newspaper.fd.df %<>% mutate(fh_status = fct_relevel(fh_status, #relevel so the legend is in order
                                             c("F", "PF", "NF")))

newspaper_country_fd_plot <- ggplot(newspaper.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, fh_scale), #reorder based on fh score
                      color=fh_status
  )) +
  labs(x = "country (listed in descending order by Freedom House score)", y = "difference in probability of supporting LGBTQs (95% CI)",
       caption = "F = Free, PF = Partly Free, NF = Not Free") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' newspaper consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray50", "gray80")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
newspaper_country_fd_plot

ggsave(filename = "figures/newspaper_country_ml_fd_plot.pdf", #Save it to directory
       plot = newspaper_country_fd_plot,
       width = 8, height = 6)

# a plot of INTERNET FD ...
# Figure A.5d
internet.fd.df %<>% subset(!is.na(fh_scale)) # remove sao tome b/c it doesn't have FH score
internet.fd.df %<>% mutate(fh_status = fct_relevel(fh_status, #relevel so the legend is in order
                                             c("F", "PF", "NF")))
internet_country_fd_plot <- ggplot(internet.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, fh_scale), #reorder based on fh score
                      color=fh_status
  )) +
  labs(x = "country (listed in descending order by Freedom House score)", y = "difference in probability of supporting LGBTQs (95% CI)",
       caption = "F = Free, PF = Partly Free, NF = Not Free") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' internet consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray50","gray80")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
internet_country_fd_plot

ggsave(filename = "figures/internet_country_ml_fd_plot.pdf", #Save it to directory
       plot = internet_country_fd_plot,
       width = 8, height = 6)

# a plot of SOCIAL MEDIA FD ...
# Figure A.5e
social_media.fd.df %<>% subset(!is.na(fh_scale)) # remove sao tome b/c it doesn't have FH score
social_media.fd.df %<>% mutate(fh_status = fct_relevel(fh_status, #relevel so the legend is in order
                                             c("F", "PF", "NF")))
social_media_country_fd_plot <- ggplot(social_media.fd.df) + 
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper, 
                      x = reorder(country, fh_scale), #reorder based on fh score
                      color=fh_status
  )) +
  labs(x = "country (listed in descending order by Freedom House score)", y = "difference in probability of supporting LGBTQs (95% CI)",
       caption = "F = Free, PF = Partly Free, NF = Not Free") +
  #ggtitle("Change in LGBTQ support when moving from \n 'none' to 'daily' social media consumption within each country") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::percent, breaks = c(-.5, -.4, -.3, -.2, -.1, 0, .1, .2, .3, .4, .5)) +
  scale_color_manual(name = "FH Status",
                     values = c("black", "gray50", "gray80")) +
  geom_hline(yintercept=c(0), linetype="dotted") + 
  theme_bw() + theme(panel.border=element_blank(), 
                     panel.grid.minor = element_blank(),
                     panel.grid.major = element_blank(),
                     plot.title = element_text(hjust=0.5)
                     #legend.position = c(.9,.5)
  )
social_media_country_fd_plot

ggsave(filename = "figures/social_media_country_ml_fd_plot.pdf", #Save it to directory
       plot = social_media_country_fd_plot,
       width = 8, height = 6)
